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ABSTRACT 


Sensor  nodes  in  a  wireless  sensor  network  (WSN)  can  establish  a  link  with  a 
UAV  by  using  beamforming  techniques  to  form  a  random  array  with  position  errors.  The 
position  errors’  effect  in  the  array  performance  is  examined  using  a  MATLAB-based 
simulation  model. 

In  order  to  spread  the  processing  and  communication  load  among  the  nodes,  two 
new  distributed  algorithms  for  beamforming  in  WSN,  based  on  the  least  squares  (LS) 
approximation  of  the  desired  array  response,  are  proposed.  The  first  is  a  distributed 
implementation  of  the  QR  decomposition,  and  the  second  is  an  iterative  method  for 
solving  the  LS  problem.  Results  indicate  that  the  processing  load  is  effectively  shared 
among  the  nodes.  Especially,  in  the  second  approach,  the  processing  load  can  be  lower 
than  that  of  the  centralized  approach,  depending  on  the  algorithm’s  convergence.  For 
both  algorithms,  the  tradeoff  for  the  ability  to  spread  the  processing  load  is  the  increased 
communication  cost,  which  could  cause  an  overall  increase  in  the  total  power 
consumption  in  the  network.  However,  the  average  power  per  participating  sensor  node  is 
still  lower  than  that  required  by  the  cluster  head  in  the  centralized  approach. 
Consequently,  the  network’s  susceptibility  to  failures  due  to  excessive  power 
consumption  is  greatly  reduced. 
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EXECUTIVE  SUMMARY 


A  wireless  sensor  network  (WSN)  consists  of  a  large  number  of  microsensors, 
each  having  limited  battery  lifetime  and  restricted  communication  and  computing 
capabilities.  Recent  advances  in  the  integrated  circuit  technology  have  allowed  the 
production  of  lightweight  and  inexpensive  sensor  nodes,  which  have  a  range  of 
capabilities,  such  as  sensing,  processing  and  communication.  WSNs  can  have  many 
applications  in  both  commercial  and  military  environments.  A  set  of  sensor  nodes,  which 
can  be  deployed  easily  and  quickly  by  an  unmanned  aerial  vehicle  (UAV),  for  example, 
can  be  used  for  monitoring  the  battlefield  environment,  sensing  for  a  wide  range  of 
targets,  especially  in  the  case  where  the  area  of  interest  is  inaccessible  or  there  is  high 
risk  of  human  loss. 

Once  deployed,  the  sensor  nodes  can  collect  the  desired  information  and  transmit 
it  to  the  UAV.  Although  a  single  sensor  node  cannot  transmit  its  data  directly  to  the  UAV 
due  to  the  limited  range  of  coverage,  several  of  them  can  coordinate  their  transmissions 
in  order  to  form  an  array  and  thus  substantially  increase  the  range  of  coverage.  Since  the 
sensors  are  randomly  deployed,  it  is  unlikely  that  they  form  topologies  that  permit  the 
formation  of  arrays  with  equally  spaced  elements.  As  a  result,  there  are  position  errors 
with  respect  to  an  array  with  equally  spaced  elements. 

A  simulation  model  was  developed  in  the  MATLAB  environment  to  analyze  the 
effect  of  these  position  errors  on  the  array  performance.  The  simulations  showed  that  the 
sidelobe  levels  in  the  array  response  increase  as  a  function  of  the  error  in  the  element 
location.  Specifically,  as  the  mean  deviation  from  the  ideal  position  was  increased,  the 
average  sidelobe  magnitudes  also  increased.  These  results  are  in  agreement  with  the 
theoretical  analysis  of  random  arrays,  found  in  the  literature.  The  degradation  of  the  array 
performance  can  be  largely  eliminated  using  a  Least  Squares  (LS)  beamformer,  which 
computes  the  weights  that  best  approximate  a  given  desired  response.  This  beamformer 
can  also  efficiently  suppress  potential  interfering  signals  coming  from  directions  other 
than  the  signal’s. 
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Since  reliability  and  robustness  in  a  sensor  network  environment  are  desired,  the 
processing  load  must  be  effectively  spread  among  the  sensor  nodes.  Centralized 
approaches  assign  the  entire  processing  load  to  a  single  node  whereas  in  a  distributed 
approach,  the  processing  tasks  are  split  into  smaller  processes,  which  are  then  allocated  to 
the  participating  sensor  nodes.  Two  fully  distributed  approaches  to  beamforming  in  WSN 
were  presented  in  this  work,  and  they  are  both  based  on  the  LS  approximation  of  the 
desired  response.  The  first  is  a  distributed  implementation  of  the  QR  decomposition,  and 
the  second  is  an  iterative  method  of  computing  the  weights  in  the  LS  sense. 

The  performance  of  the  distributed  methods  was  compared  to  the  centralized  LS 
approach  using  the  processing  and  communication  costs  as  metrics.  The  results  indicate 
that  the  processing  load  is  effectively  shared  among  the  nodes.  Especially  in  the  second 
method,  the  processing  load  is  a  function  of  the  algorithm’s  convergence  and  can  be 
lower  compared  to  that  of  the  centralized  approach,  subject  to  the  speed  of  convergence. 
For  both  algorithms,  the  tradeoff  for  the  ability  to  spread  the  processing  load  is  the 
increased  communication  cost,  which  could  cause  an  overall  increase  in  the  total  power 
consumption  in  the  network.  This  total  power,  however,  is  shared  among  the  sensor 
nodes;  therefore,  the  average  power  expended  by  a  participating  sensor  node  in  the 
distributed  implementation  is  lower  than  the  power  required  by  the  cluster  head  in  the 
centralized  approach.  Consequently,  the  network’s  susceptibility  to  failures  due  to 
excessive  power  consumption  is  greatly  reduced. 


I.  INTRODUCTION 

A.  INTRODUCTION  TO  WIRELESS  SENSOR  NETWORKS 

A  wireless  sensor  network  (WSN)  consists  of  a  large  number  of  microsensors, 
each  having  limited  battery  lifetime  and,  therefore,  restricted  communication  and 
computing  capabilities  [1],  [2]  .  Recent  advances  in  the  integrated  circuit  technology  have 
allowed  the  production  of  lightweight  and  inexpensive  sensor  nodes,  which  have  a  range 
of  capabilities,  such  as  sensing,  processing  and  communication.  If  they  are  properly 
networked  and  programmed,  these  sensor  nodes  can  cooperate  in  order  to  perform 
complex  signal  processing  functions  [3],  [4], 

The  main  issue  for  a  WSN  is  to  prolong  its  operational  lifetime  as  much  as 
possible,  taking  into  account  the  sensors’  power  consumption  requirements.  Stringent 
energy  limitations  are  also  a  crucial  factor  when  designing  signal  processing  algorithms 
for  a  WSN  [5],  Since  such  energy  restrictions  are  not  taken  into  account  by  the  signal 
processing  methods  that  are  already  used  in  applications  other  than  WSNs,  existing 
techniques  should  be  modified  in  order  to  conform  to  the  sensor  nodes’  specific 
characteristics.  Therefore,  a  major  challenge  in  recent  research  is  the  design  of  signal 
processing  and  networking  operations,  which  optimize  the  tradeoff  between  energy 
efficiency,  simplicity,  and  overall  performance,  [3]. 

Because  microsensors  are  becoming  cheaper  and  more  capable,  WSNs  will  find 
more  applications  in  both  commercial  and  military  environments  [1],  [2],  Future  tactical 
operations  will  involve  the  deployment  of  large-scale  WSNs  in  which  hundreds  or 
thousands  of  disposable  sensor  nodes  will  cooperate  in  order  to  achieve  the  mission 
objective  [3].  These  nodes  can  be  deployed  easily  and  quickly  by  an  unmanned  aerial 
vehicle  (UAV),  for  example,  as  in  Figure  1,  which  minimizes  the  risk  of  human  loss. 
Then  they  can  be  used  for  monitoring  the  battlefield  environment,  sensing  for  a  wide 
range  of  targets,  such  as  biological,  radioactive,  nuclear,  chemical  and  other  materials  [1]. 
Once  deployed,  the  sensor  nodes  can  collect  the  desired  information  and  disseminate  it  to 
a  relay  node,  such  as  a  UAV.  Furthermore,  taking  into  account  that  the  sensing 


1 


environment  may  be  harsh  and  inaccessible  for  deploying  wired  networks,  there  is 
obvious  need  for  developing  WSNs  consisting  of  small  and  disposable  sensors. 


Figure  1.  WSN  deployed  over  an  area  of  interest  and  UAV  collecting  the  desired 
information. 


B.  RELATED  WORK  IN  BEAMFORMING  AND  WIRELESS  SENSOR 

NETWORKS 

After  forming  an  ad  hoc  network  and  collecting  the  required  data  about  the  target 
of  interest,  the  sensor  nodes  must  establish  communication  with  a  UAV,  so  the  acquired 
information  can  be  transmitted  to  the  UAV.  Although  promising,  today’s  technology  still 
imposes  strict  limits  on  the  processing  and  communication  capabilities  of  the  sensor 
nodes  [6],  [7],  Single  sensor  nodes  do  not  have  sufficient  power  to  communicate  with  an 
overflying  UAV.  Since  the  UAV  may  be  required  to  fly  at  a  high  altitude  due  to  the 
hostile  nature  of  the  operative  environment,  the  objective  of  transmitting  the  collected 
information  to  the  UAV  becomes  more  difficult. 

Although  a  single  sensor  node  cannot  transmit  its  data  directly  to  the  UAV  due  to 
the  limited  batter  power,  several  of  them  can  cooperate  in  order  to  function  as  a  large 

antenna  array  and  thus  substantially  increase  the  transmission  range  and  the  data  rate. 
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This  process  of  combining  the  signal  from  different  antenna  elements  in  order  to  form  a 
single  output  of  the  sensor  array  is  known  as  beamforming.  It  has  been  proven  by  many 
studies  that  when  an  antenna  array  is  properly  configured,  it  can  improve  the  channel 
capacity  and  extend  the  range  of  coverage  [8],  [9],  It  can  also  reduce  the  multipath  fading 
and  the  bit  error  rate  (BER);  therefore,  it  results  in  more  reliable  communication  [8], 
Additionally,  beamforming  can  adaptively  steer  the  antenna  beam  towards  the  UAV,  thus 
aiming  the  radiated  energy  in  the  desired  direction  [10].  Furthermore,  the  antenna  gain  is 
proportional  to  the  number  of  the  antenna  elements,  so  the  main  beam  peak  power 
density  can  be  of  several  orders  of  magnitude  higher  than  that  of  a  single  sensor  [9], 
Another  useful  characteristic  of  an  antenna  array  is  that  it  can  be  used  in  order  to  perform 
spatio-temporal  filtering,  thus  suppressing  potential  interference  signals  coming  from 
directions  other  than  the  desired  direction  [11],  [12],  In  summary,  taking  the  above 
mentioned  advantages  into  account,  beamforming  in  WSNs  can  meet  the  objective  of 
establishing  an  efficient  communication  link  between  a  WSN  and  a  UAV. 

Several  algorithms  for  beamforming  exist  in  the  literature  [13],  [11],  [14]  and 
many  of  them  are  successfully  implemented  in  conventional  antenna  arrays.  [10], 
Nevertheless,  these  algorithms  for  the  computation  of  the  weights  for  the  array  elements 
cannot  be  directly  implemented  in  WSNs  since  there  are  significant  differences  between 
WSNs  and  conventional  arrays.  For  instance,  the  phased  arrays  used  in  RADAR  are 
installed  permanently  on  site  [15];  thus,  the  positions  of  the  antenna  elements  are  fixed. 
On  the  other  hand,  the  sensor  nodes  in  a  WSN  are  usually  randomly  deployed  and  their 
relative  positions  are  not  predetermined.  Due  to  this  random  deployment,  there  are 
position  errors,  which  cause  performance  deterioration  of  the  antenna  array,  compared  to 
that  of  an  array  of  equally  spaced  elements.  Moreover,  the  sensors  are  prone  to  frequent 
failures  due  to  limited  battery  life  or  due  to  their  vulnerability  to  environmental 
conditions.  Therefore,  the  topology  of  the  sensor  array  can  change  substantially  as  new 
nodes  are  added  or  withdrawn. 

Another  significant  problem  is  that,  in  conventional  arrays,  where  power  is  not  a 
major  issue,  the  beamforming  operation  is  performed  in  a  single  processor  [16],  All 
necessary  information  is  collected  in  a  central  processing  node,  which  is  responsible  for 
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solving  the  beamforming  problem.  However,  in  a  sensor  network  environment,  this  is 
neither  reliable  nor  desirable  since  a  single  node  would  be  assigned  this  computationally 
demanding  task.  Additionally,  such  centralized  implementations  create  a  single  point  of 
failure,  which  in  turn  creates  a  serious  system  vulnerability.  If  this  node  fails,  then  the 
beamforming  problem  has  to  be  solved  from  the  beginning.  Thus,  the  processing  load  or 
consequently  the  power  usage  should  be  effectively  distributed  across  the  sensor  network 
[17].  However,  an  optimum  set  of  participating  sensors  in  the  array  has  to  be  defined 
since  the  communication  cost  for  organizing  the  sensor  nodes  into  an  array  is  prohibitive 
after  a  certain  critical  number  of  nodes  [18]. 


C.  THESIS  OBJECTIVE 

The  objective  of  this  research  is  to  implement  several  distributed  beamforming 
algorithms,  and  to  evaluate  their  performance.  Throughout  this  work,  the  operational 
scenario  of  Figure  1  is  adopted  where  several  sensor  nodes  try  to  communicate  with  a 
UAV.  The  beamforming  process  is  not  performed  in  a  central  node  (cluster  head),  but  it 
is  split  into  smaller  processes,  which  then  can  be  allocated  to  the  sensor  nodes.  The  main 
concept  is  that  the  processing  and  communication  cost  must  be  shared  among  the  nodes, 
so  there  is  no  single  point  of  failure  and  that  the  energy  of  the  nodes  is  efficiently  used, 
thus  extending  the  WSN  lifetime.  This  work  is  focused  on  investigating  beamforming 
schemes  that  have  the  same  performance  as  well  established  centralized  approaches  yet 
offer  the  advantage  of  implementation  in  a  distributed  fashion  thus,  increasing  the 
network’s  robustness  and  overall  performance. 


D.  PROPOSED  APPROACH  TO  DISTRIBUTED  BEAMFORMING 

Starting  from  a  centralized  approach  to  the  Least  Squares  (LS)  solution  of  the 
beamforming  problem  where  the  beamformer  is  designed  in  such  way  that  the  desired 
array  performance  is  best  approximated,  two  fully  distributed  methods  are  proposed.  The 
first  one  is  a  distributed  implementation  of  the  QR  decomposition  with  Householder 
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transformations  and  derives  the  exact  solution  for  the  array  weight  vector.  The  second 
scheme  is  an  iterative  method  for  solving  the  LS  problem,  implemented  in  a  distributed 
fashion. 

In  order  to  examine  the  proposed  techniques,  a  performance  analysis  of  the 
communication  and  computational  costs  is  developed.  These  two  costs  are  closely 
connected  to  the  power  consumption  and  provide  a  reliable  test  for  the  algorithms’ 
effectiveness.  The  resulting  array  response  is  examined  and  compared  with  the  desired 
response.  The  results  from  these  two  distributed  implementations  are  encouraging  and 
indicate  that  they  can  provide  a  realistic  solution  for  the  beamforming  problem  in  sensor 
networks. 


E.  THESIS  OUTLINE 

Chapter  II  introduces  the  fundamental  concepts  of  the  antenna  arrays,  including  a 
description  of  the  uniform  linear  and  planar  array.  This  is  followed  by  an  analysis  of  the 
effects  of  position  errors  on  the  performance  of  the  antenna  array  and  a  simulation,  which 
confirms  the  theoretical  results.  Beamforming  in  wireless  sensor  networks  is  also 
presented  along  with  a  specific  operational  communication  scenario  which  uses  a  UAY. 
Finally,  the  framework  for  evaluating  the  algorithms’  performance  based  on  the  factors 
that  affect  power  consumption,  such  as  the  processing  and  communication  cost,  is 
developed. 

Chapter  III  presents  two  centralized  beamforming  approaches  and  evaluates  their 
performance.  Their  advantages  and  disadvantages  are  discussed,  and  their  ability  to 
mitigate  position  errors  is  analyzed  using  a  simulation  model  developed  in  a  MATLAB 
environment. 

In  Chapter  IV,  two  distributed  algorithms  for  beamforming  in  wireless  sensor 
networks  are  proposed.  Their  performance  is  analyzed  in  terms  of  efficiently  sharing  the 
processing  and  communication  cost  among  the  nodes,  and  the  simulation  results  are 
compared  to  those  obtained  by  the  centralized  approaches. 
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Chapter  Y  summarizes  the  significant  results  of  this  thesis  and  provides  some 
ideas  for  extending  this  work  in  the  future. 

Finally,  the  Appendix  includes  the  MATLAB  code  used  in  the  simulations. 
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II.  BEAMFORMING  IN  SENSOR  ARRAYS 


In  this  chapter,  the  main  concepts  of  antenna  arrays  are  discussed;  specifically,  the 
uniform  linear  and  planar  arrays  are  presented  as  well  as  their  array  response 
(beampattem).  The  uniform  array  is  followed  by  an  analysis  of  the  random  array  with 
position  errors.  The  effects  of  the  antenna  element  position  errors  in  the  main  lobe  and 
the  sidelobe  power  gain  are  presented  under  various  assumptions  about  the  statistical 
characteristics  of  the  error.  Beamformer  implementations  are  described,  particularly  the 
narrowband  beamformer.  Lastly,  there  is  a  discussion  about  communication  and 
computational  cost  in  sensor  networks  as  a  function  of  power  consumption,  which 
indicates  the  need  for  distributed  algorithms  in  beamforming. 


A.  UNIFORM  LINEAR  AND  PLANAR  ARRAYS 

The  fundamental  concepts  of  the  uniform  linear  and  planar  arrays  are  presented  in 
this  section  along  with  the  basic  formulation  of  beamforming  which  will  be  used 
throughout  this  work. 

1.  One-Dimensional  Array 

In  Figure  2,  a  linear  array  is  depicted  with  M  identical,  equally  spaced  elements. 
The  spacing  between  consecutive  array  elements  or  sensor  nodes  in  the  case  of  a  sensor 
array  is  assumed  equal  to  a  half  wavelength,  i.e.,  d  =  A/2,  where  the  elements  are 
isotropic,  meaning  their  beampattem  is  omni-directional.  Furthermore,  it  is  assumed  that 
the  array  is  located  far  enough  from  the  signal  source  (e.g.,  a  UAV);  thus,  it  is  considered 
to  lie  in  the  source’s  far  field.  The  array  axis  is  assumed  to  be  in  the  x-axis.  The  plane 
wave  s(t )  arrives  at  the  array  at  an  angle  0a  with  respect  to  the  x-axis  (array  axis),  and 

for  this  reason,  the  m  +  Ith  element  senses  the  signal  earlier  than  the  mth  element.  If  xm  is 

the  distance  between  the  m,h  node  and  the  reference  node,  which  is  located  at  the  origin 
of  the  coordinate  system,  then  the  signal  s{t )  arrives  at  the  m‘h  element  earlier  by  tm 
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seconds  with  respect  to  the  reference  element.  The  time  difference,  which  depends  on  the 
arrival  angle  0  or  Angle  of  Arrival  (AO A),  and  the  element’s  distance  from  the 

reference  point  is  given  by  [19]: 


/  x  xm  cos  6 l 

K)-- — ' 


where  c  is  the  speed  of  light. 

z-axis 


X2 

< - ► 

X3 

Xm 

Figure  2.  An  Mxl  linear  array  of  equally  spaced  isotropic  elements. 

Each  array  element  is  weighted  by  a  complex  weight  wm,  for  m  =  0,1,..., M, 
which  multiplies  the  incoming  signal.  Adding  all  the  elements’  weighted  inputs  gives  the 
spatial  response  of  the  array  or  array  factor  F(6a )  for  any  arbitrary  angle  Ba : 

M 

=  (2) 

m= 1 

where  wm  =  1  and  I m  and  e'0>'"e‘"  are  the  magnitude  and  the  phase  of  the 
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complex  weights,  respectively.  Using  the  wavenumber  /?  =  2k  I X  and  the  expression  for 

the  time  difference  tm,  the  weights  can  be  written  as  wm  =  Jmej/3XmC0Sda<>  and  the  array 
factor  as 

M 

(3) 

m= 1 

The  weights  wm  are  carefully  selected  in  order  to  give  the  maximum  value  of  the 
array  response  F(6a )  at  the  desired  direction  0d  and  to  suppress  potential  interference 
signals  arriving  from  other  directions.  Indeed  at  0a  =  0  ,  the  array  response  reaches  its 
maximum  value 

M  M 

F(0aJ  =  fjIme-j^cos^ej^c°^  =X/ffl  =M  if  Im  =  l,Vm  (4) 

m= 1  m= 1 

The  set  of  weights  wm  form  the  weight  vector 

w  =  [wl  w2  ...  wmf 

while  the  steering  vector  or  direction  vector  is  defined  as 

)  =  [1  eJfiX2COS0a  e  iPXi  C0S°"  eifixhta*9a  ]r 

and  incorporates  the  location  information  of  the  array.  Therefore,  the  array  response  can 
be  expressed  as 

F{Oa)  =  wH  d{0a).  (5) 


Note  that  the  main  concept  of  beamforming  is  the  use  of  the  weights  w  in  order  to  point 
the  array  beam  towards  any  desired  direction.  So,  if  the  desired  transmission  direction  is 
0  ,  then  the  beamformer  should  set  its  weights  to  be 


w 


=  7 


jj3xm  cos  6a 


(6) 


The  beampattem  of  a  10x1  uniform  linear  array  is  shown  in  Figure  3  where  the 
normalized  power  gain  G  is  defined  as  [9] 
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(7) 


GWa)  = 


\mf 

max|F(^)|2 


and  it  is  plotted  as  a  function  of  the  direction  0a  (in  degrees).  The  array  elements  are 
identical  and  isotropic,  and  the  spacing  has  a  fixed  value  of  A I  2  and  beam  pointing 
angle  0  =  0° .  The  maximum  sidelobe  is  equal  to  -13  dB  and  the  Half  Power  Beamwidth 

(HPBW)  is  about  10.2°. 


Figure  3.  Normalized  Power  Gain  (dB)  -  Beampattem  of  a  M  =  10  element  array  with 
isotropic  elements,  fixed  spacing  A  /  2  and  0  =  0° . 

In  general,  the  sidelobe  level  decreases  with  increasing  number  of  elements  M 
and  approaches  the  value  of -13.3  dB.  The  HPBW  in  any  plane  containing  the  array  axis, 
is  given  [9]  by 

0MB  =  0.866 


Md 
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where  Md  □  A ,  i.e.,  it  is  valid  for  long  arrays.  Therefore  the  HPBW  is  a  function  of 
M  and  decreases  as  the  number  of  array  elements  increase. 


2.  Two-dimensional  (planar)  Array 

The  previous  discussion  of  a  linear  array  can  be  easily  expanded  to  a  planar  array, 
and  similar  expressions  can  be  derived  for  the  array  response  and  the  power  gain.  In 
Figure  4,  a  uniformly  spaced  planar  array  is  depicted  with  MxN  identical  and  isotropic 
elements.  The  spacing  between  the  array  elements  is  assumed  equal  to  half  wavelength, 
d  =  A/2,  in  both  directions  while  it  is  assumed  again  that  the  array  is  located  in  the 
source’s  far  field.  The  plane  wave  s(t)  arrives  at  the  array  at  polar  angle  0O  with  respect 

to  the  z-axis  and  an  azimuth  angle  (j)o  with  respect  to  the  x-axis;  thus,  the  ( m,n)th  element 
receives  the  signal  earlier  by  tmn  seconds  compared  to  the  reference  element  at  the  origin. 
This  time  difference  tmn  depends  on  the  angles  0() ,  (f>0  and  the  element’s  position 
( xmn->ymn )  in  the  array,  and  is  given  by  [19] 


Knn  (<?0.A) 


nn  Sin  P  C0S  <K  +  /Ajnn  sin  «£■  <K 

C 


(8) 


Adding  all  the  elements’  weighted  inputs  gives  the  array  response  of  the  planar 
array  F{0,  <f>)  for  any  arbitrary  choice  of  angles  0  and  (f)  : 


N  M 


*  „jP(xmn  sin^cos^+jOTW  sin^sin^) 


w  e 

mn 


n= 1  m= 1 


(9) 


where 


w  —  J  V,„  sin  C°S'A|  sin  ^  sin  (j,H  ) 

vmn  ±  mn’' 


(10) 


are  the  complex  weights  for  each  element  ( m ,  n) .  In  matrix  form,  the  above  expression 
can  be  written  as 


F(0,<f>)  =  wH  d(0,</>) 


(11) 


where  w  is  an  MN  x  1  weight  vector  and 
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1 

gjP(x 12  sin cos (/)+yn  sin#sin^) 
gj/3(x\T,  sin^cos^+j13  sin#sin^) 


ejfi(x\M  sin6>cos^+j1M  sin  #  sin  ^) 


(12) 


eJPixMN  sin  6>  cos  (f)+yMN  sintfsin^) 


is  an  MN  x  1  steering  vector. 

The  weights  wmn  are  again  selected  in  order  to  give  the  maximum  value  of  the 
array  response  F{0,(p )  in  the  desired  direction (#0  (j>{]) .  At  6  =  6()  and  (f)  =  (f>{) ,  the  array 
response  reaches  its  maximum  value 

F(0o,<po)  =  MN  if  Imn=l,Vm,n  (13) 

where  in  (13)  it  is  assumed  that  the  excitation  of  elements  (weighting)  is  uniform. 
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z-axis 


a)  M  xN  planar  array  of  equally  spaced  isotropic  elements 


b)  Measurement  of  angles  for  direction  of  signal  arrival  (d0,</>0) 


Figure  4.  Signal  wavefront  arriving  from  direction  at  a  M  x  N  planar  array. 
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In  Figure  5,  the  three-dimensional  beampattem  of  a  10x10  uniform  planar  array 
is  shown,  where  the  normalized  power  gain  G  in  dB  is  plotted  as  a  function  of  the  polar 
angle  6  and  the  azimuth  angle  (f) .  A  cross  section  of  this  3-D  beampattem  is  plotted  in 
Figure  6  where  the  azimuth  angle  is  constant  at  <f>{]  =45°  (direction  of  arrival),  and  the 
beampattem  varies  with  the  polar  angle  6 . 
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Figure  5.  Normalized  Power  Gain  (dB):  3-D  Beampattem  of  a  10x10  planar  array  with 
isotropic  elements,  equally  spaced  by  A/2  and  direction  of  signal 
arrival  (#0  =  30°,  <p0  =  45°) . 
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Figure  6.  Normalized  Power  Gain  (dB):  Beampattem  of  a  10x10  planar  array  with 
isotropic  elements,  equally  spaced  by  A/2  and  direction  of  signal 
arrival  (0{)  =  0°,  </)()  =  45°) .  Azimuth  angle  is  fixed  at  </>0  =  45° . 


The  maximum  sidelobe  is  equal  to  -26  dB,  and  the  HPBW  is  about  10.2°.  In 
general,  the  sidelobe  levels  decrease  with  an  increasing  number  of  elements  MN  and 
approaches  the  value  of  -26.6  dB.  Similarly,  the  HPBW  decreases  continuously  as  the 
number  of  elements  increase. 


B.  RANDOM  ARRAY:  POSITION  ERRORS 

In  a  random  deployment  of  a  sensor  array  where  the  sensor  nodes  are  dropped 
randomly  over  an  area  of  interest,  it  would  be  unrealistic  to  expect  formations  of 
perfectly  spaced  planar  arrays.  An  illustrating  example  comes  from  the  fact  that  for  an 
operating  frequency  fc  =  1  GHz,  the  ideal  distance  between  the  sensor  nodes  is 

A/2  =  15  cm,  and  even  a  small  displacement  of  3  cm  yields  a  position  error  of  20%. 
Therefore,  the  performance  of  a  sensor  array  should  be  studied  using  the  theoretical 
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analysis  of  random  arrays.  In  general,  the  beampattem  of  a  randomly  deployed  sensor 
array  will  be  affected  by  position  errors,  amplitude  and  phase  errors,  quantization  errors 
and  failures  of  the  nodes,  which  cause  a  change  in  the  array  topology.  Throughout  this 
work,  the  main  emphasis  will  be  given  to  the  effect  of  position  errors  and  solutions  to  this 
problem. 


1.  Position  and  Phase  Errors 


There  are  several  references  [19],  [20],  [21],  [22]  in  the  literature  that  deal  with 
position  errors  in  random  arrays.  The  random  antenna  elements’  misplacement  causes 
phase  errors  and  mismatches,  which  yield  degraded  performance  for  the  array.  In  [19],  an 
analysis  of  the  radiation  pattern  of  a  random  array  with  both  amplitude  and  phase  errors 
is  presented.  For  an  M  xN  two-dimensional  array,  assuming  that  there  are  no  amplitude 
errors,  i.e.,  the  weights  wmn  have  the  same  magnitude  I  =Imn,  the  expected  increase  As  in 
the  sidelobe  level  with  respect  to  the  main  lobe  is  given  by  [19] 


A„  = 


(14) 


where  the  phase  errors  follow  a  Gaussian  distribution  with  zero  mean  and  variance  a\  . 
For  example,  in  a  5x5  planar  array,  A^  is  equal  to  0.00518  for  oA(p  =20°  and  equal  to 
0.025  for  oA(p  =  40° .  Therefore,  doubling  of  the  phase  error  will  cause  an  increase  of 

almost  6  dB  in  the  sidelobe  level  with  respect  to  the  main  lobe.  However,  it  can  be  seen 
from  (14)  that  the  effect  of  phase  errors  can  be  mitigated  by  increasing  the  number  of 
array  elements  in  the  array.  Indeed,  As  can  be  decreased  by  a  factor  of  two  if  the  number 
of  antenna  elements  is  doubled. 

The  fractional  loss  in  the  main  lobe  gain  due  to  phase  errors  is  given  by  [19] 

—  =  e-<  (15) 

G„ 


where  G  and  G0  are  the  main  lobe  power  gains  with  and  without  the  presence  of  phase 
errors,  respectively.  Since  for  a  3  dB  reduction  in  the  main  lobe  level,  a  phase  error 
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standard  deviation  of  a ^  =  47°  is  needed,  it  can  be  concluded  that  the  main  lobe  is  not 
significantly  affected  by  random  position  errors. 

2.  Random  Array  Implementation 

In  general,  the  effects  of  misplacement  errors  depend  strongly  on  the  assumptions 
about  the  random  characteristics  of  the  position  deviations  from  the  uniform  array.  In  the 
previous  section,  the  results  were  derived  based  on  the  assumption  of  Gaussian  phase 
errors.  However,  there  are  several  references  in  the  literature  which  consider  different 
deployments  of  sensor  arrays  and  consequently  different  assumptions  about  the  random 
distribution  of  the  phase  errors.  One  such  example  is  included  in  [22]  where  the  location 
of  each  node  in  the  sensor  array  is  chosen  randomly  by  a  uniform  distribution  within  a 
disk. 

Throughout  this  work,  the  position  errors  will  be  modeled  as  a  deviation  from  the 
ideal  position  of  a  uniform  array.  The  position  errors  will  be  modeled  as  a  uniformly 
distributed  random  variable  with  a  minimum  value  of  0  and  a  maximum  value  of 
ax  1/2  where  a  is  the  maximum  percentage  error;  therefore,  the  mean  error  will  be 
a  x  1  /  4 .  It  is  also  assumed  that  there  is  displacement  in  both  x  and  y  directions  as 
indicated  in  Figure  7.  The  array  axis  for  randomly  placed  elements  is  the  best  line  fit  of 
the  sensor  nodes’  positions  and  it  is  assumed  to  be  also  the  x-axis  of  the  coordinates 
system.  In  Figure  7  the  best  line  fit  of  the  nodes  location  is  found  and  it  is  assumed  to  be 
the  array  axis.  Then  the  deviations  are  defined  using  this  array  axis  as  a  reference.  The 
beampattem  is  computed  in  a  plane  which  is  perpendicular  to  the  x-y  plane  and  at  <j)(]  (in 

this  case  (f>()  =  45°)  with  respect  to  the  x-axis. 

In  Figure  8,  the  effect  of  position  errors  in  a  10x1  linear  antenna  array  topology 
is  depicted.  The  sidelobe  levels  have  been  increased  and  consequently  the  array 
performance  has  been  deteriorated.  The  mean  beampattem  of  a  10x1  linear  array  with 
mean  position  deviation  equal  to  0.3/1/ 2  in  both  x  and  y  directions  is  compared  to  the 
beampattem  of  the  ideal  uniform  10x1  array  for  polar  angle  6{]  =  30°  and  azimuth  angle 

<f>i}  =  45°.  The  array  axis  is  in  the  x-direction  and  the  beampattem  is  in  a  plane 
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perpendicular  to  the  x-y  plane  and  at  45°  with  respect  to  x-axis  ($,  =45°).  For  the 

random  case,  the  mean  beampattem  is  obtained  after  averaging  the  beampattems  for  50 
repetitions  of  randomly  generated  array  topologies. 
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Figure  7.  Ideal  and  approximately  linear  arrays  with  randomly  deployed  elements  and 
position  errors. 
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Figure  8.  Beampattem  of  a  10x1  linear  array  with  position  errors,  compared  to  a  10x1 
array  with  equally  spaced  isotropic  elements.  Polar  angle  0O  =  30° ,  azimuth  angle 

(j)0  =  45°  and  mean  deviation  30%  of  A  /  2  in  both  x-y  directions. 

Next,  the  relationship  between  the  position  errors  and  the  increase  in  the  average 
sidelobe  levels  is  plotted  in  Figure  9.  The  deviation  from  the  ideal  position  is  modeled  as 
a  uniform  random  variable  in  the  range  of  0  to  A/2  in  both  x  and  y  directions,  i.e., 
mean  deviation  is  0.5A/2.  For  each  value  of  mean  position  error,  the  increase  for  the 
first  sidelobe  (red  line)  and  the  largest  of  all  the  sidelobes  (blue  line)  is  calculated  by 
averaging  50  simulation  runs.  It  is  obvious  that  the  sidelobe  power  gain  increased 
significantly  as  the  mean  deviation  increased;  for  example,  if  the  mean  deviation  is  about 
0.4T/ 2,  then  the  maximum  sidelobe  increase  is  almost  5.3  dB  while  the  first  sidelobe 
increase  is  about  4.2  dB.  Thus,  the  strongest  sidelobe  is  lower  from  the  main  lobe  by  7.7 
dB  only  while  the  first  sidelobe  differs  from  the  main  lobe  by  8.8  dB. 
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Figure  9.  Effect  of  position  errors  in  a  random  array.  Average  sidelobe  level  (dB)  for 
the  first  (red  line)  and  the  largest  (blue  line)  sidelobe  as  a  function  of  the  mean 
deviation  of  the  actual  element  positions  from  the  ideal  positions  in  linearlOxl 
array  with  equally  spaced  isotropic  elements. 

C.  BEAMFORMER  IMPLEMENTATIONS 

In  the  previous  sections,  it  was  shown  that  the  beampattem  of  an  antenna  array  is 
determined  by  the  direction  of  the  incoming  signal,  the  number  of  the  array 
elements MxN ,  and  the  array  topology,  which  includes  the  position  errors  and  the  set  of 
weights  wmn .  The  objective  of  a  beamformer  is  to  preferentially  receive  a  signal  from  a 

specific  direction  or  to  preferentially  transmit  a  signal  in  that  direction.  It  is  also  usually 
desirable  to  suppress  interference  signals,  which  come  from  other  directions.  Therefore, 
the  beamforming  operation  consists  of  adjusting  the  weights  wmn  in  such  way  that  the 

main  lobe  is  steered  towards  the  desired  signal’s  AO  A.  Such  a  beamformer  for  M 
elements  in  a  linear  array  is  depicted  in  Figure  10,  and  it  is  typically  used  for  narrowband 
signals. 
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Figure  10. 


The  output  of  the  beamformer  is  given  by 

y{t)  =  wH  x{t)  (16) 

where  w  is  the  weight  vector  and  x(t)  is  the  signal  vector 

x(t)  =  [x,  (t)  x2(t)  ...  (17) 

A  conventional  beamformer  can  be  described  [8]  as  a  delay-and-sum  beamformer 
with  all  weights  having  the  same  magnitude.  The  phases  are  selected  to  steer  the  array 
beampattem  towards  a  desired  direction  (#0,^0).  However,  there  are  many  types  of 

beamformers,  which  can  be  classified  as  either  data  independent  or  statistically  optimum, 
depending  on  how  the  weights  are  chosen  [11]. 

In  a  data  independent  beamformer,  the  weights  are  chosen  so  as  to  create  a 
specified  desired  response  for  all  signal  and  interference  cases.  The  array  data  (the  signal 
vector  x(t) )  are  either  not  known  or  not  taken  into  account  for  the  beamforming  design. 
If  the  desired  response  is  Fd  (6,  (f>)  in  order  to  receive  a  signal  from  a  certain  direction  and 
cancel  out  interferences  from  other  directions,  then  the  weight  vector  w  is  chosen  in  such 
a  way  that  the  actual  response  F(0,</>)  approximates  the  desired  one.  In  the  following 
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chapters,  the  beamforming  operation  will  be  based  on  data  independent  techniques, 
which  try  to  create  an  approximation  of  the  desired  response  Fd  {6,  </>) . 

The  second  class  of  beamformers  [11]  contains  the  statistically  optimum  ones.  In 
this  case,  the  weights  are  chosen  based  on  the  statistics  of  the  array  data,  and  the  goal  is 
to  optimize  the  array  response  using  several  criteria.  Generally,  a  statistically  optimum 
beamformer  tries  to  cancel  out  the  interfering  signals  by  placing  nulls  in  their  incoming 
directions  in  order  to  maximize  the  Signal  to  Noise  Ratio  (SNR)  at  the  output  of  the 
beamformer.  Such  a  beamformer  is  the  Multiple  Sidelobe  Canceller  (MSC),  which  needs 
auxiliary  channels  free  of  the  desired  signal  but  is  very  simple  in  its  implementation. 
Another  optimum  beamformer  requires  knowledge  of  both  the  desired  signal  and  noise 
covariance  matrices  Rs  and  Rn ,  respectively,  and  has  the  advantage  that  it  maximizes  the 

SNR,  but  it  is  not  easy  to  implement.  The  Minimum  Variance  Distortionless  Response 
(MVDR)  beamformer  computes  the  weights  given  a  set  of  constraints  and  provides  a 
very  good  performance,  but  it  is  also  computationally  intensive. 

There  are  also  adaptive  algorithms  for  beamforming,  which  compensate  for  the 
fact  that  the  signal’s  statistics  are  usually  not  known  or  may  vary  with  time.  Such 
beamformers  for  the  weight  determination  are  based  on  the  well-defined  and  popular 
LMS  and  RLS  algorithms  and  also  on  numerous  variations  of  them.  A  more  thorough 
analysis  of  these  adaptive  algorithms  and  of  the  previously  described  statistically 
optimum  beamformers  can  be  found  in,  [1 1],  [12]  and  [14] . 


D.  WIRELESS  SENSOR  NETWORK  ARRAYS 

This  section  presents  a  discussion  of  sensor  node  deployment  schemes  and  the 
concepts  of  the  processing  and  communication  costs  which  will  be  used  for  the 
performance  evaluation  of  the  beamforming  algorithms. 

1.  Deployment  of  Sensor  Nodes 

As  mentioned  in  Chapter  I,  recent  developments  in  the  MEMS  technology  have 


enabled  the  construction  of  small,  cheap,  multifunctional  sensors  with  signal  processing 
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and  communication  capabilities.  These  sensor  nodes  can  be  deployed  over  an  area  of 
interest  and  can  be  used  in  order  to  collect  process  and  transmit  information.  However, 
there  are  still  many  challenges  for  the  efficient  implementation  of  a  sensor  network,  such 
as  efficient  power  consumption,  controllable  deployment  of  the  nodes,  source 
localization,  self  organization,  and  others. 

In  this  work,  the  specific  application  scenario  which  is  examined  is  summarized 
in  Figure  1 1 .  A  set  of  sensor  nodes  are  deployed  in  the  battlefield,  and  they  acquire  the 
desired  information,  which  may  be  any  type  of  signal  (acoustic,  video,  etc.).  The  UAV, 
flying  over  the  sensor  field,  establishes  a  connection  with  the  network  in  order  to  obtain 
the  collected  data.  Due  to  the  limited  sensor  capabilities,  a  single  node  can  not  transmit 
its  data  directly  to  the  UAV  since  it  does  not  have  the  required  transmission  range. 
However,  the  nodes  can  be  organized  into  a  large  antenna  array  where  each  sensor  plays 
the  role  of  an  antenna  element  and  implements  a  distributed  beamformer. 


Figure  11.  A  sensor  network  deployment  used  for  information  transfer  to  a  UAV  (After 
Ref.  [34]). 
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In  the  previous  section,  it  was  shown  that  the  best  performance  for  an  antenna 
array  is  achieved  when  the  elements  are  located  on  a  rectangular  grid  with  an 
interelement  distance  equal  to  A  /  2 .  Nevertheless,  it  is  obvious  that  this  ideal  deployment 
can  rarely  be  achieved  in  a  real  operational  environment.  The  nodes  may  be  dropped  by 
the  UAV  or  deployed  by  a  ground  force  and  consequently  the  sensor  array  topology 
cannot  be  the  desired  one.  This  randomness  imposes  position  errors,  which  effect  the 
beampattern  of  the  sensor  array  as  analyzed  in  the  previous  section. 

Much  work  has  been  done  in  analyzing  the  dependence  of  the  random  array 
response  on  the  statistical  characteristics  of  the  array  topology  [19],  [20],  [21], [22], 
Furthermore,  much  emphasis  has  been  placed  on  algorithms  that  allow  a  node  to  define 
its  position  with  respect  to  a  reference  node  within  the  sensor  network  [23].  In  the  next 
chapter,  several  beamforming  algorithms  that  use  location  information  will  be  described. 
Assuming  that  the  coordinates  of  a  node  in  a  local  coordinate  system  are  known,  there  are 
many  schemes  for  successfully  mitigating  the  effect  of  the  position  errors  and  choosing  a 
suitable  set  of  weights  for  the  beamformer.  However,  there  are  limitations  in  the 
performance  of  these  algorithms  when  the  position  errors  are  large  enough. 

Another  practical  approach  could  be  to  search  and  locate  a  suitable  subset  of 
nodes  and  then  form  an  antenna  array.  In  [24]  and  [25],  a  central  node  (cluster  head)  tries 
to  find  a  set  of  nodes  whose  topology  is  close  to  a  uniform  linear  or  planar  array 
according  to  some  geometric  criteria.  Using  only  the  distance  between  the  nodes  as 
known  information,  the  proposed  algorithm  tries  to  find  the  optimum  subset  of  sensors 
with  minimum  mean  position  error.  Figure  12  illustrates  this  concept,  where  a  set  of  five 
sensors  yielding  the  best  approximation  to  a  5x1  uniform  linear  array  is  determined. 
After  finding  this  optimum  set  of  nodes,  several  beamforming  schemes  can  be  applied  in 
order  to  find  the  weight  coefficients.  Therefore,  the  combination  of  finding  a  suitable  set 
of  nodes  with  an  appropriate  beamforming  technique  can  provide  sufficient  performance. 
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Figure  12.  Finding  five  nodes  approximating  a  5x1  linear  array  with  equally  spaced 
elements  (After  Ref.  [24]). 


2.  Communication  and  Computational  Cost 

An  important  issue  in  the  sensor  networks  is  the  efficient  energy  consumption  by 
the  nodes.  Low-power  hardware  components  and  low-duty  cycle  operation  techniques 
must  be  applied  in  order  to  achieve  the  required  ultra-low-power  operation.  Energy  is 
consumed  while  processing  or  transmitting  data.  Therefore,  the  beamforming  algorithms 
must  be  evaluated  in  terms  of  realistic  power  consumption.  In  order  to  achieve  this  goal, 
it  is  obvious  that  there  is  a  need  for  schemes  implemented  in  a  distributed  manner;  thus, 
the  computational  effort  and  consequently  the  energy  needed  is  shared  among  the 
sensors.  Moreover,  these  techniques  should  minimize  the  communication  since  wireless 
transmission  consumes  considerable  amount  of  power  during  a  node’s  operation. 

The  need  for  distributed  algorithms  in  order  to  minimize  the  communication 
power  consumption  is  discussed  in  [26].  A  radio  transmitting  1  kb  of  data  over  a  distance 
of  100  m,  with  an  operating  frequency  of  1  GHz  using  BPSK  modulation  having  an  error 
probability  of  10'6  and  fourth-power  distance  loss  with  Rayleigh  fading,  requires 
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approximately  3  Joules  of  energy.  The  same  amount  of  energy  can  perform  300  million 
instructions  for  a  100  MIPS/watt  general  processor  [26].  This  results  in  a  ratio  of  30,000 
processing  instructions  per  transmitted  bit  with  equal  energy  consumption.  Other 
practical  implementations  [26],  [27]  have  yielded  ratios  from  200  to  3000.  This  ratio  of 
communication  cost  to  computational  cost  depends  largely  on  the  sensor  characteristics 
(transmission  range,  complexity  of  the  instruction  in  number  of  bits,  etc.)  and  can  be 
increased  in  the  presence  of  noise  where  retransmissions  will  be  needed. 

In  general,  the  relationship  of  the  communication  and  computational  cost  with  the 
power  consumption  depends  on  the  technical  specifications  of  the  sensors,  the 
applications,  and  other  factors  that  are  difficult  to  exactly  predict,  such  as  the  presence  of 
noise.  However,  a  general  framework  can  be  formed  using  the  definitions  in  Table  1. 


Symbol 

Definition 

P* 

Mean  power  per  transmitted  bit  (power  consumed  to  transmit  a  number  of 

bits  divided  by  this  number) 

Ntb 

Number  of  transmitted  bits 

Pi 

Mean  power  per  instruction  in  the  sensor’s  processor  (power  consumed  to 

perform  a  number  of  instructions,  divided  by  this  number) 

N, 

Number  of  instructions 

Pc 

Total  transmission  power  or  communication  load 

Pp 

Total  processing  power  or  computational  load 

p 

Total  power  for  a  specific  application 

Table  1 .  Quantities  used  for  defining  the  communication  and  computational  cost. 
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The  transmission  power  Pc  is  given  by 


Pc=NtbxPtb 

and  the  processing  power  is  given  by 

Pp=NlxPl 

while  the  total  power  is  the  sum  of  these  two  factors 

P  =  Pc+Pp- 


(18) 


(19) 


(20) 


The  parameter  77  is  defined  as  the  ratio  of  the  power  per  transmitted  bit  to  the  power  per 

instruction,  and  as  mentioned  before,  it  can  vary  from  200  to  3000,  indicating  that  one 
transmitted  bit  requires  much  more  energy  than  one  performed  instruction  by  the  sensor’s 
processor. 


n<„  = 


P» 


(21) 


Throughout  this  work  the  communication  cost  is  focused  on  the  implementation 
of  the  algorithms.  Therefore,  it  depends  only  on  the  data  elements  that  need  to  be 
transmitted  for  the  implementation  of  the  different  beamforming  schemes.  These  data 
elements  must  be  encapsulated  in  data  packets,  thus  there  is  a  communication  overhead 
due  to  these  packets.  Furthermore,  the  organization  of  the  sensor  nodes  in  a  cluster  needs 
also  a  number  of  control  packets.  This  networking  cost,  which  depends  on  various  factors 
as  the  number  of  sensors,  the  noise  and  the  protocols  used,  will  not  be  taken  into  account 
throughout  this  work  and  will  be  left  for  future  analysis. 

Using  the  above  definitions,  the  total  power  consumption  for  the  implementation 
can  be  calculated  and  the  different  beamforming  methods  can  be  compared  using  a 
common  framework. 


27 


E.  SUMMARY 


In  this  chapter,  the  basics  of  beamforming  in  antenna  arrays  were  discussed; 
specifically,  the  beampattems  of  the  uniform  linear  and  planar  arrays  were  presented. 
This  was  followed  by  an  analysis  of  the  effect  of  the  position  errors  in  random  arrays  on 
the  array  performance  as  measured  by  an  increase  in  the  sidelobe  levels.  Beamformer 
implementations  were  discussed,  including  a  brief  reference  to  various  techniques  that 
have  been  reported  in  the  literature.  The  specific  operational  scenario  for  communication 
between  sensor  networks  and  UAV  was  also  briefly  presented.  Finally,  the 
communication  and  computational  cost  as  functions  of  the  power  consumption  were 
introduced  as  metrics  for  the  evaluation  of  beamforming  algorithms  in  sensor  networks. 
The  strict  requirements  for  low  power  consumption  by  the  sensor  nodes  create  the  need 
for  distributed  beamformimg  algorithms  (presented  in  Chapter  IV)  compared  to 
centralized  algorithms  (presented  in  Chapter  III). 
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III.  CENTRALIZED  IMPLEMENTATION  OF  A  BEAMFORMER 


This  chapter  is  focused  on  centralized  implementation  of  the  beamforming 
operation  for  linear  arrays.  In  the  centralized  approach,  it  is  assumed  that  all  information 
needed  to  determine  the  weight  vector  is  available  in  a  specific  node,  which  can  be  the 
cluster  head  in  a  sensor  network.  This  node  will  collect  all  the  data,  such  as  the  steering 
vector  d_(6,  <f>) ,  which  depends  on  the  array  topology,  the  direction  of  the  desired  signal 
(#„,$,)  and  the  direction  of  potential  interferences,  and  will  calculate  the  weight  vector 
to  provide  the  desired  array  response  at  the  output  of  the  beamformer. 

Two  different  implementations  are  presented  and  their  advantages  and 
disadvantages,  such  as  simplicity  in  implementation,  overall  performance  and 
computational  demands,  are  discussed.  Finally,  an  analysis  of  the  computational  and 
communication  cost  examines  the  feasibility  of  these  implementations  in  sensor  networks 
where  power  efficiency  is  a  major  issue. 

A.  PHASE  MATCH  OF  THE  STEERING  VECTOR 

This  is  the  simplest  implementation  of  a  beamformer  and  can  be  considered  an 
extension  of  the  conventional  beamformer  for  the  case  of  random  position  errors.  In  this 
method,  as  in  the  case  of  a  uniform  array  (see  Chapter  II),  the  weights  of  the  elements  w 
are  chosen  with  the  goal  to  match  the  steering  vector  d(d0)  and  create  the  main  lobe  with 
maximum  gain  towards  the  AOA  O0 .  If  the  steering  vector  is 

=  eJPx2sine0  eifixs si"  eiPxn  sin^ojr  (22) 

then  the  weights  are  w  =  d_(d0) .  Each  element  of  the  weight  vector  has  unit  magnitude 

and  the  same  phase  as  the  corresponding  element  of  the  steering  vector.  Thus,  the  array 
response 

F(d)  =  wH  d(d)  (23) 

has  its  maximum  F(0O)  at  the  AOA  6>0 . 
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Due  to  position  errors  in  the  sensor  array,  the  distances  x1,x2...,xn  from  the 
reference  node  are  not  multiples  of  the  ideal  distance  XU.  However,  the  weight  vector 
can  still  be  selected  to  match  the  steering  vector  in  the  desired  direction  of  6() .  The 

implementation  of  this  approach  in  a  sensor  network  is  very  simple  and  assumes  the 
following: 

a)  Each  node  can  calculate  its  position  from  a  reference  node  accurately. 

b)  The  relative  positions  of  the  nodes  are  disseminated  to  a  central  processing 
node,  a  cluster  head  if  a  cluster  hierarchical  architecture  is  established. 

c)  The  desired  steering  direction  of  the  incoming  signal  is  known  to  the  central 

node. 

d)  Errors  due  to  noise  during  the  communication  among  the  nodes  are  not  taken 
into  account. 

The  cluster  head  (central  node)  collects  all  the  necessary  information  and 
calculates  the  weight  vector.  Then  it  sends  to  each  node  its  corresponding  weight,  which 
will  be  used  for  the  beamforming  and  the  steering  of  the  sensor  array. 

The  array  pattern  of  this  simple  beamformer  is  shown  in  Figure  13  for  a  10x1 
linear  array  with  uniform  position  errors  up  to  0.4X  /  2 .  The  signal’s  direction  of  arrival  is 
known  ( 60  =30°,  =45°).  The  figure  shows  the  mean  beampattem  averaged  over  50 

simulation  runs. 

Some  conclusions  can  be  drawn  for  this  beamforming  technique.  First,  it  is  simple 
in  implementation  since  the  central  processing  node  needs  only  to  receive  the  relative 
position  from  each  node  one  at  a  time  and  then  send  back  the  calculated  weights  to  each 
node.  Second,  the  main  lobe  of  the  beampattem  has  the  maximum  value  in  the  direction 
of  arrival  (#„,$,) ,  and  it  is  exactly  equal  to  the  maximum  value  of  the  (ideal)  uniform 

array.  However,  the  beampattem  presents  notable  deviations  for  angles  other  than  the 

AOA.  In  this  specific  case,  the  sidelobe  power  level  increases  significantly  for  angles 

lower  than  0°.  In  general,  the  performance  of  the  array  approaches  that  of  the  uniform 

array  at  angles  around  the  AOA  but  deteriorates  for  other  elevation  angles. 
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Figure  13.  Mean  beampattem  of  a  random  10x1  uniform  linear  array  with  position 
errors,  averaged  over  50  simulation  runs,  compared  to  a  uniform  10x1  array  of 
equally  spaced  elements.  The  AOA  is  ( 6()  =30°,  <f>()  =45°).  The  position  errors 
are  uniformly  distributed  between  0  and  and  0.4/1  /  2  in  both  x  and  y  directions. 
Beamforming  is  implemented  by  matching  the  steering  vector. 

B.  BEAMPATTERN  APPROXIMATION  IN  THE  LEAST  SQUARE  SENSE 

The  beamformer  presented  in  the  previous  section  computes  the  weights  of  the 
antenna  elements  by  trying  to  achieve  the  best  performance  in  the  direction  of  the  signal 
arrival  and  does  not  set  any  requirement  for  the  other  directions.  Although  it  behaves  well 
at  angles  around  the  desired  direction,  it  suffers  from  high  sidelobes  in  other  directions. 
Therefore,  it  cannot  be  used  when  there  are  interference  signals  coming  from  directions 
other  than  the  signal’s  AOA  or  if  the  desired  signal  comes  from  a  certain  range  of 
directions.  If  there  is  a  strong  source  of  interference,  the  desired  array  response  should  be 
zero  in  that  direction,  in  order  to  cancel  the  interfering  signal.  As  a  result,  there  is  a  need 
for  a  beamforming  design,  which  calculates  the  weight  coefficients  w  in  such  way  that 
the  resulting  response  approximates  the  desired  response  over  a  range  of  directions. 
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1.  Least  Squares  Problem  Formulation 

As  mentioned  before,  the  array  response  for  a  specific  angle  6  is  given  by 

F(6)  =  wH  d{6).  (24) 

If  the  desired  response  F_d{6 )  is  defined  over  m  number  of  angles,  then  the  actual  array 
response  F_(0)  of  n  antenna  elements  for  these  m  angles  is 

F(0)H  =  whD(9 )  (25) 

where 

=  F(02)  -  F(6m)]T  (26) 

is  a  fiix  1  vector  which  contains  the  array  response  for  each  one  of  the  angles  0i ,  with 
i  =  1, . . .  m .  The  weight  vector  w  is  an  n  x  1  column  vector  and  D(6)  is  an  n  x  m  steering 
matrix  containing  the  steering  vectors  d_(0:)  for  each  one  of  the  angles  9t : 


II 

Q 

m)  - 

m,)\ 

(27) 

ejP*\  sin  <9, 

gj P x\  sin  02 

ejpx\  sin  6m 

D{6)  = 

ejpx2  sin6». 

ejpx2  sin 

gjpx2  sin6»w 

(28) 

5‘ 

ejPxn  sin<92 

gjpXn  sin  0m 

Note  that  for  simplicity  the  array  response  is  defined  over  a  set  of  arrival  angles 
while  the  azimuth  angle  <f>{)  is  fixed.  If  the  desired  response  is  defined  over  a 

combination  of  m  elevation  angles  6v...,0m  and  q  azimuth  angles  </>^...,</>q,  then  the 
array  response  F{6 ,  (f>)  is  a  m  /  q  matrix  and  the  steering  matrix  D(6,  (f>)  is  an  nxmxq 
3-D  matrix.  This  adds  unnecessary  complexity  to  the  formulation  of  the  problem  and  is  to 
be  avoided.  In  the  following  analysis,  the  array  response  will  be  defined  for  a  known, 
fixed  azimuth  angle  and  a  set  of  m  arrival  angles 
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If  the  desired  response  is  F_d  ( 0 ) ,  the  goal  is  to  select  the  weights  w  such  that  the 
actual  response  Fid)  generated  by  this  set  of  weights  approximates  the  desired  response. 
This  can  be  done  by  trying  to  minimize  the  L2  -norm  of  the  difference  between  the 
desired  and  the  actual  response.  In  other  words,  the  weight  coefficients  w  are  chosen  in 
order  to  minimize  the  mean  squared  error  between  F_d  ( 6 )  and  Fid) 

s  =  rmn\F(S)-FA0)\2.  (29) 


Therefore  the  problem  of  the  antenna  weight  calculation  is  solved  using  a  classic  least 
squares  (LS)  procedure.  Taking  the  transpose  of  (25),  the  array  response  can  be  written  as 

F(6)  =  D(0)hw.  (30) 


For  m  number  of  angles  or  approximation  points  and  n  number  of  sensors,  with  m>n  , 
the  problem  becomes  an  overdetermined  LS  problem 


s 


=  min 


DiO)Hw-Fd 


(31) 


where  again  F(d)  is  a  m  x  1  vector  and  D{6)H  is  a  m  x  n  matrix.  Provided  that  the  n  x  n 
matrix  D{6)D{6)H  is  invertible  (i.e.,  D{6)H  is  full  rank),  then  the  solution  to  the  LS 
problem  of  (30)  is  given  by 

™  =  D+(0)Fd(0)  (32) 

where  D+  (0)  is  the  pseudo  inverse  of  D{0)H  defined  as 

D+(0)  =  (D(0)D(0)HylD(0).  (33) 

The  weight  coefficients  that  have  been  calculated  using  the  LS  approach  have  both 
different  amplitudes  and  phases.  In  the  array  with  equally  spaced  elements,  the 
amplitude  Im  was  the  same  for  all  elements  ( I m  =  1  ,Vm )  and  only  the  phase  was  different. 

However,  in  the  LS  approach  the  weights  have  also  different  amplitudes.  So  in  the  LS 
approaches  the  sensors  modify  both  amplitudes  and  phases  in  order  to  approximate  the 
desired  response. 
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The  design  of  an  array  using  the  above  decribed  LS  approximation  of  the  desired 
array  response  is  related  to  the  design  of  arrays  using  the  Fourier  analysis.  In  this 
approach,  the  desired  array  function,  which  is  required  to  create  the  desired  array 
response,  can  be  expanded  in  a  Fourier  series.  Truncating  this  Fourier  series,  the  result  is 
an  array  with  a  finite  number  of  elements.  Such  an  array  is  optimum  in  the  sense  that  no 
other  array  with  the  same  number  of  elements  can  approximate  the  desired  array  function 
with  lower  mean  squared  error. 

2.  Implementation  and  Performance  Analysis 

The  significant  advantage  of  the  above  LS  approach  is  that  the  desired  response 
can  be  closely  approximated  over  a  set  of  angles,  thus  providing  the  ability  to  cancel  out 
unwanted  signals  and  amplifying  only  the  desired  one.  Results  of  such  an  implementation 
of  the  LS  scheme  for  a  10x1  random  array  can  be  seen  in  Figure  14,  where  the  desired 
response  Fd  (0)  is  the  array  response  of  the  (ideal)  uniform  10x1  linear  array.  For  the 
approximation  of  the  desired  beampattem  180  points  (angles  from  -90°  to  90°)  were 
used.  The  mean  beampattern  is  averaged  over  50  simulation  runs,  and  there  is  virtually  no 
difference  between  the  desired  (red  line)  and  the  actual  beampattem  (black  line)  as 
generated  by  the  weights  computed  by  minimizing  the  LS  criterion.  The  actual  response 
can  be  also  compared  with  the  response  derived  by  the  technique  of  the  phase  match 
approach  from  the  previous  section  (cyan  line). 
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Figure  14.  Mean  beampattem  of  a  random  10x1  linear  array  with  position  errors  (blue) 
and  mean  beampattem  calculated  using  the  “phase  match”  (cyan)  and  the  LS 
approach  (red),  averaged  over  50  simulation  runs,  compared  to  a  uniform  10x1 
array  (black).  The  AOA  is  (  90  =  30° ,  </>0  =  45°  ).  The  position  errors  are  uniformly 
distributed  between  0  and  and  0.5/1/ 2  in  both  x  and  y  directions 


The  performance  is  almost  the  same  for  both  phase  match  and  LS  schemes  around 
the  desired  AOA,  but  the  LS  approach  is  obviously  better  in  any  other  direction.  The 
phase  match  technique  is  more  vulnerable  to  interference  signals,  especially  to  those 
coming  from  0  <  0° ,  while  the  LS  approach  creates  a  beampattem  almost  identical  to  the 
desired  response  of  a  uniform  array.  Indeed,  the  determination  of  the  weights  using  the 
LS  approximation  technique  achieves  elimination  of  the  effect  of  the  position  errors  (blue 
line);  the  position  errors  are  uniformly  distributed  between  0  and  0.5/1/ 2  in  both  x 
and  v  directions. 

The  LS  approach  is  clearly  a  centralized  approach.  The  implementation  in  a 
sensor  network  makes  similar  assumptions  as  in  the  phase  match  approach  including  an 
additional  assumption:  the  desired  array  response  defined  over  a  set  of  angles  0v...,9m  is 
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assumed  to  be  known  to  the  cluster  head.  This  response  may  include  the  direction  (#„,$,) 
of  the  incoming  signal  and/or  the  directions  of  interfering  sources. 

The  cluster  head  collects  all  the  necessary  information  and  specifically  the 
relative  positions  xl,...,xn  of  the  nodes  and  the  desired  response  F d (6) .  The  steering 
matrix  D(0)  is  completely  defined  by  the  positions  and  the  desired  angles  for  the 
beampattem  approximation.  After  collecting  all  necessary  data,  the  cluster  head 
constructs  the  steering  matrix  using  the  position  data  and  the  set  of  angles  and  calculates 
the  weight  vector  by  solving  the  LS  problem  of  (30).  Each  node  receives  the 
corresponding  coefficient  from  the  cluster  head,  and  the  sensor  array  is  ready  to  transmit 
or  receive  towards  the  desired  signal’s  direction  and  simultaneously  suppress  potential 
interfering  sources. 

The  LS  approach’s  performance  is  satisfactory  in  eliminating  the  effect  of  the 
position  errors  as  observed  in  Figure  14.  The  assumptions  mentioned  before,  such  as  that 
information  about  the  relative  position  of  the  nodes  and  the  desired  response  (incoming 
direction  of  signal  and  interferences),  are  reasonable  in  a  sensor  network  environment. 
There  are  algorithms  which  allow  a  sensor  node  to  calculate  its  relative  position  from  a 
cluster  head  using  localization  techniques  [23],  [26].  Thus,  the  LS  scheme  could  be 
implemented  in  a  cluster-based  sensor  array  where  all  information  is  collected  and 
processed  in  the  cluster  head. 

Nevertheless,  this  centralized  approach  suffers  from  the  inherent  problems  of  any 
scheme  where  the  processing  effort  is  not  distributed  among  the  nodes  but  performed  by 
a  single  node  (cluster  head).  Since  the  main  effort  is  undertaken  by  one  sensor  only,  a 
single  point  of  failure  is  created  in  the  sensor  array.  If  this  specific  node  fails,  then  the  LS 
problem  has  to  be  solved  from  the  beginning,  i.e.  a  new  central  sensor  must  collect  all  the 
information,  calculate  the  weight  coefficients  and  disseminate  them  to  the  other  nodes. 
This  is  a  waste  of  valuable  processing  and  communication  resources  which  are  very 
restricted  in  the  sensor  nodes.  Considering  that  failure  of  sensor  elements  in  a  sensor 
array  is  something  common  due  to  their  low  battery  life  and  low-cost  construction,  this 
centralized  approach  lacks  robustness. 
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Furthermore,  since  new  nodes  may  offer  to  participate  in  the  sensor  array  or  other 
nodes  may  decide  to  switch  into  sleep  mode  due  to  limited  power,  the  array  topology  will 
change  frequently  and  dramatically.  This  means  that,  for  any  modification  in  the  sensor 
array,  the  LS  problem  must  be  solved  from  the  start  in  order  to  determine  the  new  set  of 
weight  coefficients.  For  this  reason,  the  centralized  approach  does  not  offer  the  required 
scalability,  which  is  desired  in  a  sensor  network  environment. 

Finally,  the  processing  load  of  solving  the  LS  problem  increases  exponentially  as 
the  number  of  sensors  n  and  the  number  of  approximation  points  (angles)  m  increase.  If 
this  processing  effort  is  performed  by  a  single  sensor,  then  its  power  will  be  consumed 
rapidly  and  the  sensor  will  fail. 

The  solution  is  a  distributed  implementation  of  the  LS  approach,  and  such  two 
schemes  will  be  presented  in  Chapter  IV.  First,  the  single  point  of  failure  limitation  can 
be  overcome  by  distributing  the  processing  load  among  many  nodes.  The  total  processing 
effort  in  a  distributed  approach  may  be  higher  compared  to  the  centralized  approach,  but 
the  tradeoff  is  the  increased  robustness  of  the  system.  Furthermore,  the  LS  approach  must 
be  implemented  in  such  way  that  changes  in  the  array  topology  do  not  require  the 
solution  of  the  problem  from  scratch  but  only  small  modifications  to  the  already 
calculated  weight  coefficients.  This  can  be  achieved  by  an  algorithm  implemented  in  a 
distributed  fashion,  which  will  offer  considerable  scalability  to  the  problem.  Since  the  LS 
approach  performs  very  well  in  approximating  the  desired  pattern,  a  distributed 
implementation  of  it  will  first  retain  this  feature  and  second,  it  will  offer  the  desired 
scalability  and  robustness;  thus,  it  can  provide  a  reliable  solution  for  beamforming  in 
sensor  networks. 

In  the  next  section,  a  discussion  about  the  processing  cost  as  a  function  of  the 
number  of  approximation  points  is  presented.  As  mentioned  before,  the  processing  load 
increases  dramatically  with  the  number  of  approximation  points;  therefore,  the  minimum 
number  of  points  should  be  used  in  order  to  reduce  the  computational  effort. 
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c. 


APPLICATION  OF  LEAST  SQUARES  BEAMFORMING  APPROACH  TO 
SENSOR  NETWORKS 


In  the  previous  section  the  LS  approach  was  presented  for  the  approximation  of  a 
desired  response  defined  over  a  set  of  direction  angles  or  approximation  points.  As 
mentioned  before,  the  processing  load  of  the  LS  approach  depends  largely  on  the  number 
of  nodes  and  the  number  of  approximation  points.  For  example,  if  n  is  the  number  of 
sensors  and  m  the  number  of  angles,  the  solution  of  the  LS  problem  using  the  normal 

equations  needs  (m  +  ry^)n2  flops  [28],  By  using  fewer  approximation  points,  the 

processing  load  can  be  reduced,  which  is  crucial  for  the  sensor  nodes  with  limited 
processing  power.  However,  the  solution  will  not  be  the  same  as  the  unmodified  case, 
and  the  approximation  will  degrade  as  the  number  of  points  is  decreased.  Therefore,  it  is 
worth  examining  the  “quality”  of  the  desired  response  approximation  as  a  function  of  the 
number  of  approximation  points. 

In  order  to  illustrate  this  relationship,  a  series  of  50  simulations  with 
corresponding  random  arrays  was  performed.  In  each  run,  the  desired  response  of  a 
uniform  linear  10x1  array  (ideal)  was  approximated  using  different  number  of 
approximation  points  ranging  from  10  to  90.  For  each  set  of  approximation  points,  the 
corresponding  solution  for  the  weight  coefficients  was  derived  using  the  LS  approach. 
The  reference  approximated  response  and  the  reference  weight  vector  was  based  on  180 
approximation  points  (one  for  each  degree  from  -90°  to  90° ).  The  solution  of  the  LS 
approach  for  each  set  of  approximation  points  is  compared  to  the  reference  solution 
derived  from  180  points. 

Two  error  metrics  are  used  to  evaluate  the  performance  of  each  set  of 
approximation  points.  The  first  metric  ew  is  the  mean  L2  -norm  of  the  difference  between 
the  reference  weights  and  the  examined  weights 
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where  wm  are  the  weights  derived  from  the  180-point  approximation  (reference),  w.  are 
the  weights  derived  by  an  i-  point  approximation,  for  /  =  10,...,  90,  and  S  =  50  is  the 
number  of  simulation  runs. 


The  second  metric  eF  is  the  mean  L2  -norm  of  the  difference  between  the 

reference  array  response  derived  from  1 80  approximation  points  and  the  array  response 
for  fewer  approximation  points  (10  to  90) 


SF 
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where  Fm) {6)  is  the  approximated  array  response  using  180  points  (reference)  and  Ft{0) 
is  the  response  derived  by  an  i  -  point  approximation,  for  i  =  10, . . . ,  90 . 

The  test  topology  is  a  10x1  random  array  where  the  antenna  elements  have 
position  errors,  modeled  as  uniform  random  variables  between  0  and  0.4/1  /  2 .  The  results 
are  shown  in  Figures  15  and  16  for  the  two  errors  ew  and  eF ,  respectively. 


Figure  15.  Error  metric  ew  of  weights  as  a  function  of  the  number  of  approximation 
points  m  . 
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Figure  16.  Error  metric  eF  of  approximated  array  responses  as  a  function  of  the  number 
of  approximation  points  m  . 

In  Figure  15,  error  metric  ew  decreases  as  the  number  of  approximation  points 

increases,  which  is  reasonable.  The  significant  result  comes  from  the  fact  that  the 
difference  between  the  reference  weights  and  the  weights  computed  for  fewer 
approximation  points  is  very  small,  even  for  15  approximation  points.  This  means  that 
using  only  15  approximation  points,  the  weight  vector  is  very  close  to  the  reference 
weight  vector  (180  approximation  points).  Therefore,  there  is  no  need  to  use  180  points  to 
approximate  the  desired  response  since  this  can  be  done  by  using  only  15  points,  which 
implies  that  the  processing  cost  can  be  significantly  reduced. 

Similar  comments  can  be  made  about  the  results  in  Figure  16  for  the  error 
between  the  reference  array  response  derived  by  180  approximation  points  and  the 
responses  derived  by  fewer  points.  If  the  desired  response  is  known  and  needs  to  be 
approximated  using  a  LS  approach  for  the  determination  of  the  weights,  the  number  of 
approximation  points  can  be  substantially  lowered  without  significant  degradation  in  the 
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performance  of  the  antenna  array.  Thus,  the  processing  cost  can  be  reduced,  which  is 
desirable  in  a  sensor  network  environment  with  restrictions  on  the  power  consumption. 


D.  SUMMARY 

In  this  chapter,  two  centralized  approaches  for  the  computation  of  the  weights 
were  presented.  The  first  technique  is  simple  and  easy  to  implement,  but  the  performance 
is  not  satisfactory  in  the  presence  of  interference  signals.  The  second  approach  is  based 
on  the  LS  approximation  of  a  known  desired  response  by  selecting  the  weights  in  order  to 
minimize  the  error  between  the  actual  and  the  desired  response.  Although  having 
satisfactory  performance,  the  LS  scheme  still  lacks  robustness  and  scalability  since  it  is  a 
centralized  approach.  Finally,  the  performance  of  the  approximation  of  the  desired 
response  as  a  function  of  the  number  of  approximation  points  was  examined.  The  next 
chapter  presents  two  algorithms  that  implement  the  LS  approach  in  a  distributed  fashion. 
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IV.  DISTRIBUTED  ALGORITHMS  FOR  BEAMFORMING 


A  Least-Squares  (LS)  based  centralized  approach  for  solving  the  beamforming 
problem  was  described  in  the  previous  chapter.  Nevertheless,  in  a  sensor  network 
environment  where  the  processing  and  communication  load  must  be  shared  among  the 
nodes,  there  is  a  need  for  distributed  beamforming  algorithms. 

Two  distributed  implementations  of  the  LS  beamforming  problem  are  proposed  in 
this  chapter.  The  first  is  based  on  the  well  known  QR  factorization  using  Householder 
transformations,  where  each  node  performs  a  part  of  the  QR  decomposition  process.  The 
second  approach  is  a  distributed  iterative  solution  of  the  LS  problem,  which  converges 
quickly  to  the  actual  weight  coefficients.  Both  algorithms  efficiently  distribute  the 
processing  load  among  the  nodes;  however,  the  tradeoff  consists  of  an  increase  in  the 
communication  cost. 


A.  DISTRIBUTED  QR  FACTORIZATION  WITH  HOUSEHOLDER 
TRANSFORMATIONS 

1.  Householder  Transformations 

A  Householder  transformation  (or  Householder  reflection)  is  a  transformation  of 
reflecting  a  vector  about  some  known  plane  [28],  [29],  Given  an  arbitrary  vector 

x  =  [xj  x2  •••  xm]reDm  and  a  unit  vector  e,  =  [l  0  •••  0 ^  &Um ,  the 

mxm Householder  matrix  is  defined  as  the  matrix  that  transforms  x  to  a  vector  parallel 
to  ex 


H  x  =  ae  j 


where  a  is  a  scalar. 

A  Householder  matrix  can  be  formed  by  [29] 

H  =  I-  —  vv 

P 


(36) 


(37) 
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where  /?  is  a  scalar,  /  is  an  identity  matrix  and  y  is  an  m  x  I  vector.  From  the  above 
definition,  the  Householder  matrix  H  is  completely  determined  by  vector  y  and  scalar 
/? ,  so  there  is  no  need  to  store  all  m2  elements  of  H,  but  only  y  and  fJ ,  which  are  only 
m  + 1  elements. 

The  significant  characteristic  of  the  Householder  matrix  is  that  it  is  an  orthogonal 
matrix;  thus,  it  has  the  property 

HT=Hl=H  .  (38) 

Orthogonal  matrices  or  orthogonal  transformations  can  be  used  to  obtain  a  QR 
factorization  of  a  matrix  A ,  and  this  in  turn  can  be  used  to  solve  a  linear  system  Ax  =  b, 
as  described  in  the  following  section. 

The  above  defined  Householder  matrix  zeros  out  all  the  elements  of  an  mxl 
vector  x  except  the  first  one,  but  one  can  construct  Householder  matrices  that  zero  out 

only  the  last  m-k  components  of  a  vector  x  [29].  Let  x 1 1  and  xl2>  define 

T(1)  =[xi  x2  •••  -v,  f  (39) 

*(2)=K  xm  •••  xmf  (40) 

while  1^  and  / ^  denote  {k - 1) x {k - 1)  and  ( m-k  +  l)x(m-k  +  \ )  identity  matrices, 
respectively.  From  (32),  an  {m  -  k  + 1)  x  (m  -  k  + 1)  Householder  matrix  H(k2)  can  be 
constructed  as  [29] 

wfl=rj>-i-viviT  («) 

Pk 

such  that 


Hl2)x(2)  = 


x 


.(2) 


Zi 


(42) 


By  defining  the  mxm  matrix  Hk  as 
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(43) 


Hk  = 


7(1) 

0 


H\ 


(2) 


the  following  results 


Hkx  = 
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h — H 

1 _ 

X(1) 
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**  Xr  1 

(  m  \ 

Vx2 
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•  0 
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<N 

_ i 
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(N 

<N 

_ 1 
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k-l 

V  i=k  ) 

(44) 


The  Householder  matrix  Hk  defined  in  (43)  is  an  orthogonal  matrix  that  can  be 
written  as 


Hk  =I-  —  vkvTk 

Pt 


(45) 


and  acts  like  an  identity  matrix  on  the  first  k- 1  coordinates  of  any  vector  yen'"  and 
transforms  the  rest  into  a  unit  vector.  Moreover,  it  is  not  necessary  to  store  the  entire 
matrix  Hk ,  but  it  is  enough  to  store  the  (m-k  + 1)  x  1  vector  vk  and  the  scalar  j3k . 


2.  QR  Decomposition 

The  QR  decomposition  is  used  in  order  to  factor  a  matrix  A  into  a  product  of  two 
matrices  Q  and  R  [28],  [29] 


A  =  QR 


(46) 


If  A  is  an  m  x  n  matrix,  then  Q  is  an  m  x  m  orthogonal  matrix  and  R  is  an  m  x  n  upper 
triangular  matrix.  The  QR  decomposition  is  used  to  solve  a  linear  system  Ax  =  b,  and  it 
can  be  computed  by  applying  a  series  of  Householder  transformations  to  matrix  A. 

By  defining 


hi 

an 

an  ■ 

••  <hn 

'21 

a22 

a23  ■ 

••  a2  n 

hi 

a32 

a33  • 

■■  a3  n 

(47) 
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and  using  the  procedure  described  in  the  previous  section,  an  mxm  Householder  matrix 
Hx  can  be  found,  which  when  applied  to  the  first  column  of  A  will  give  a  multiple  of  ex . 

Thus,  the  result  of  multiplying  Hx  with  A  will  be 


l~a(1) 

u\\ 

a(l) 
u  12 

aB  • 

'•  a?n 

0 

^  • 

HXA  = 

0 

a33  ' 

••  a(l) 

U3n 

0 

am2 

a(l)  ■ 

um3 

••  a(1) 

mn  J 

(48) 


where  the  superscript  denotes  that  the  matrix  element  is  affected  by  the  Hx 
transformation.  Note  that  Hx  is  fully  described  by  vector  vx  and  scalar  /?,  as  mentioned 
before.  Then,  a  second  Householder  transformation  H2  that  will  zero  out  the  last  n-  2 
elements  in  the  second  column  of  HXA  can  be  again  constructed.  Thus,  H2HXA  will  be 
of  the  form 


r  a(1) 
W11 

a(1) 

w12 

a(1)  • 

u13 

••  a(l) 

U\n 

0 

a{2) 

u22 

a(2)  • 

u23 

••  a(2) 

U2n 

H2HxA  = 
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0 

42)  • 

••  a(2) 

U3n 

0 

0 

til  • 

••  a(2) 

mn 

(49) 


where  again  the  superscript  indicates  the  last  transformation  that  affected  the 
corresponding  element.  For  example,  the  first  column  and  the  first  row  are  not  affected 
by  the  H2  transformation.  From  (41),  H2  is  of  the  form 


H2  = 


Ix  0 
0  H™ 


(50) 


where  Ix  is  a  lxl  identity  matrix,  and  H{2)  is  an  (m  - 1)  x  (m  - 1)  Householder  matrix. 


Since  H(2)  is  fully  determined  by  an  (m  -1)  x  1  vector  v_2  and  a  scalar  ,  only 


(2) 


m  elements  are  needed  to  be  stored  for  H 


(2) 
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Continuing  to  apply  the  Householder  transformation  in  this  way  will  result  in  an 
upper  triangular  matrix  R  ,  i.e., 


HnHn_l-H2HlA  =  R 

where  R  is  of  the  form 


a(l) 
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(51) 


(52) 


The  mxn  matrix  R  is  composed  of  two  matrices,  a  nxn  upper  triangular  Rx  and  a 
( m-n)xn  zero  matrix: 


R  = 


0 


(53) 


Let  Qt  ,  an  mxm  matrix,  be  the  product  of  the  Householder  transformations 


QT=HnH n_l-  H2Hx 


(54) 


where  Qx  is  an  nxm  matrix  containing  the  first  n  rows  of  QT ,  and  Q2  is  an 

(m-n)xm  matrix  containing  the  last  in  -  n  rows.  Since  QT  A  =  R  and  Q  is  orthogonal, 
it  follows  that 


A  =  QR  =  [Q{ 


Q2] 


R, 


=  QA 


(55) 
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The  above  derived  QR  decomposition  can  be  used  in  order  to  find  the  solution  to 
the  least  squares  problem  Ax  =  b,  where  A  is  an  mxn  matrix  and  x  and  b  are  n x  1  and 
mxl  vectors,  respectively.  The  solution  x  satisfies  the  normal  equations  [29]: 

AT  Ax  =  AT  b  (56) 


and  by  using  the  QR  decomposition,  they  can  be  written  in  the  form 


By  setting 


c  =  Qrb  = 


Qi  k 
02 


Tb 


(57) 


(58) 


and  taking  into  account  that  QXQX  =  I  and  Rj  is  nonsingular,  (57)  simplifies  to 


Rxx  =  cx  .  (59) 

The  system  of  (59)  can  be  solved  using  back  substitution  [29]  and  the  solution 

x  =  Rl~1c1  (60) 


is  the  unique  solution  to  the  least  squares  problem. 

In  summary,  if  A  is  an  mxn  matrix,  the  least  squares  problem  can  be  solved 
using  the  QR  decomposition  by  Householder  factorizations  as  follows: 

a)  Compute  R  =  HnHn_x  •  •  •  H2HXA  and  c  =  HnHn_x  ■  ■  ■  H2HX  b  . 


b) 


Partition  R  and  c  into  a  block  form  R  = 

and  c  = 

_o_ 

_<2__ 

where  Rx  and  cx 


each  have  n  rows. 

c)  Finally,  solve  Rxx_  =  cx  using  back  substitution. 
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3.  Proposed  Algorithm  Description 

In  Chapter  III,  a  method  for  determining  the  weights  that  provide  the  best 
approximation  of  a  desired  response  in  the  LS  sense  was  described.  The  LS  problem  of 
(30)  can  be  solved  using  the  QR  factorization  presented  in  the  previous  section.  The 
cluster  head  collects  all  the  positions  from  the  sensor  nodes  and  constructs  the  steering 
matrix  D(0)H  .  Then  it  computes  the  QR  factorization  using  Householder  transformations 
and  finds  the  solution  of  the  weight  vector  w  for  a  given  desired  response  F_d  {6) .  The 
procedure  is  straightforward,  but  it  lacks  robustness  since  a  single  node  bears  the  entire 
processing  load,  which  increases  the  possibility  of  failure,  causing  the  problem  to  be 
solved  again  from  scratch.  Therefore,  a  distributed  algorithm  for  solving  the  LS  problem 
with  the  QR  decomposition  is  desirable.  The  processing  load  is  shared  by  the  nodes,  but 
there  is  also  a  tradeoff  of  increased  communication  load. 

For  the  implementation  of  the  algorithm,  the  following  reasonable  assumptions 
are  made: 

a)  Each  node  can  calculate  its  position  from  a  reference  node  accurately. 

b)  Each  node  can  broadcast  information  to  other  nodes  in  the  cluster. 

c)  The  desired  array  response  for  a  set  of  m  directions  0v...,0m  of  the  incoming 
signal  is  known. 

d)  Errors  due  to  the  noise  during  the  communication  among  the  nodes  are  not 
taken  into  account. 


The  algorithm  exploits  the  specific  nature  of  the  matrix  D(0)"  to  be  decomposed 
by  the  QR  factorization.  By  taking  the  transpose  of  (28),  the  steering  matrix  is 


Dh{0)  = 


ejP*\  sin#! 
ejP*\  sin  #2 

gj Px\  sin  0rn 


gj Px2  sin  0\ 
g  j Px2  sin  #2 

ejpx2  sin  0m 


e.)Pxn  sin#, 
ejpxn  sin  #2 

e  JPxn  sin  dm 


(61) 
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where  n  is  the  number  of  sensor  nodes  and  m  the  number  of  angles  or  approximation 
points.  Obviously,  the  first  column  is  constructed  by  the  position  of  the  first  node  xx  and 

the  set  of  angles  0x,...,0m.  Similarly,  the  elements  of  the  second  column  depend  on  the 
second  node  position  x2  and  the  set  of  m  angles.  In  general,  each  column  i  of  DH  (0) 
depends  only  on  the  position  x;  of  the  corresponding  node  and  the  desired  angles. 

From  (61),  it  is  obvious  that  since  one  node  knows  its  position  from  a  reference 
node  and  the  set  of  m  angles,  it  can  construct  its  corresponding  column  of  D"  (0) 
without  exchanging  any  information,  such  as  their  position,  with  the  other  nodes.  Thus, 
the  matrix  DH  (0)  can  be  stored  in  a  distributed  way  among  the  sensor  nodes. 

From  the  analysis  of  the  previous  section  comes  the  observation  that  the  first 
Householder  transformation  Hx  is  needed  to  zero  out  all  the  elements  of  the  first  column 

except  for  the  first  element.  Since  for  the  construction  of  matrix  Hx  only  the  data  values 

from  the  first  column  are  used,  the  conclusion  drawn  is  that  it  can  be  computed  by  the 
first  node  only.  Therefore,  the  first  node,  which  initializes  the  QR  decomposition, 
calculates  the  Hx  matrix  and  multiplies  its  own  first  column  by  this  matrix.  Then  it 

broadcasts  Hx  to  all  other  nodes,  which  in  turn  transform  their  corresponding  columns 
with  Hx .  There  is  no  need  to  transmit  the  entire  matrix  Hx  but  only  the  vector  y,  and  the 
scalar  J3X.  After  this  first  step  of  Hx  computation,  transmission  and  application  is 
completed,  the  distributed  steering  matrix  will  have  the  form  of  (48): 
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where  are  the  matrix  elements  transformed  by  Hx . 

The  second  step  is  similar  to  the  first  and  is  performed  by  the  second  node.  The 

second  node  computes  the  H2  transformation,  which  depends  only  on  the  second 
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column,  and  then  H2  (vector  v2  and  scalar  /?2 )  is  broadcast  to  each  node  which  in  turn 

applies  the  new  Householder  transformation  to  its  column.  The  resulting  distributed 
matrix  is 
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(63) 


From  the  structure  of  D(H2j(0) ,  the  following  observations  are  made.  The  first 
node  receives  the  H2  transformation,  but  it  does  not  need  to  apply  it  to  its  column  since 
H2  does  not  affect  the  first  column.  Additionally,  each  node  does  not  need  to  apply  the 
H2  to  the  first  element  of  its  column  since  the  first  row  is  not  affected  by  the  H2 

transformation.  Similarly,  the  following  Householder  transformations  are  applied  only 
when  needed,  thus  avoiding  redundant  calculations. 

The  procedure  goes  on  until  all  Householder  transformations  are 

computed,  broadcast  and  applied  by  the  nodes.  The  distributed  matrix  will  result  in  an 
upper  triangular  form 
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The  QR  decomposition  of  the  DH  (< 9 )  matrix  is  achieved  in  a  fully  distributed  way,  and 
the  processing  load  is  shared  among  the  nodes.  There  is  no  additional  processing  effort 
since  the  Householder  transformations  are  applied  only  to  the  affected  elements,  as 
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described  before.  There  is  no  difference  between  the  centralized  and  the  distributed  QR 
decomposition  in  terms  of  the  processing  load.  The  significant  advantage  is  that  instead 
of  one  node  bearing  all  the  computational  burden,  each  node  in  the  array  takes  its  turn  in 
completing  the  QR  decomposition.  However,  there  is  some  additional  communication 
load  since  the  //  matrices  need  to  be  broadcast.  This  is  the  tradeoff  for  relieving  the 
cluster  head  from  the  heavy  processing  load.  This  first  phase  of  the  algorithm  is  depicted 
in  Figure  17.  In  the  1st  step  the  1st  node  calculates  the  matrix  Hx  and  broadcasts  it  to  all 

the  other  sensors.  In  the  2nd  step  the  2nd  node  calculate  the  matrix  H2  and  the  procedure 
continues  until  the  last  node  which  computes  the  last  matrix  //3 

Ist  step 


2nd  step 

Hi 

© 

3rd  step 

Hi.  H2 

Figure  17.  First  phase  of  the  algorithm  for  distributed  QR  decomposition  by  the 

sensor  nodes.  Bolded  and  underlined//;  above  the  nodes  denote  the  Householder 

transformations  that  are  already  stored  in  the  sensor.  Simple  //  denote  the  just 
computed  and  broadcast  Householder  transformation. 

The  next  phase  in  solving  the  LS  problem  is  to  update  the  desired  response  using 
a  series  of  Householder  transformations  to  obtain  vector  c : 
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(65) 


which  is  then  partitioned  to  c  = 


according  to  (58). 


The  third  phase  includes  the  solution  of  the  system 


Rxw  =  Cj  = 


(66) 


by  back  substitution  where  Rx  defined  by  (52)  and  (53)  is  an  upper  triangular  n  x  n 
matrix  containing  the  first  n  rows  of  D''n)(6) : 
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The  nth  node  needs  to  solve  the  last  equation  of  the  system 


a^w  =c 

nn  n  n 


(68) 


and  immediately  computes  its  weight  wn  for  beamforming.  The  (n  -  \)th  node  needs  to 
solve  the  ( n  -  l),h  equation  in  order  to  find  its  weight  : 

a(n~y>  ,  w  ,+a('n~})w  =  c  . .  (69) 

n-l,n-l  n- 1  n-l,n  n  n- 1  V  / 
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The  data  element  of  matrix  Rx  in  (67)  is  constructed  by  the  (n-  l)th  node  by 

applying  all  transformations  Hx  to Hn  ,  to  the  initial  element  an_ln_x  =  eiPXn-'sme"-'  since 
this  data  element  belongs  to  its  own  (n  - \)th column.  However,  it  needs  the  value , 
which  is  in  the  nth  column;  thus,  it  depends  on  the  position  xn  of  the  nth  node.  It  also 
needs  to  receive  the  weight  wn  and  compute  the  value  cn_x .  The  nth  node  broadcasts  its 
position  xn  and  its  weight  wn  to  all  nodes.  Now,  the  (n  -  \)th  node  can  solve  (69)  and 
calculate  the  weight  wn[ . 


Similarly,  the  ( n  -  \)th  node  broadcasts  its  position  xn[  and  its  weight  wn_t  to  all 
the  other  nodes,  which  need  this  information  to  solve  their  own  equation.  For  example, 
the  ( n  -  2)th  node  has  to  solve  the  equation 


a("2)  w  7+a("72).w  . +a("72)w  =c  7, 

n-2,n-2  n-2  n-2,n~\  n- 1  n-2,n  n  n-2  ? 


(70) 


so  it  needs  the  positions  xn  and  xn_x  in  order  to  construct  the  and  a(nn_~2)n  elements, 

respectively,  by  applying  all  the  Householder  transformations  up  to  Hn_ 2 .  It  also  needs 
the  weights  wn_x  and  >v(i  from  the  previous  nodes. 

The  final  result  is  that  every  node  obtains  its  own  weight,  and  the  sensor  network 
is  now  ready  to  cooperate  in  order  to  form  an  antenna  array.  The  third  phase  of  the 
algorithm  is  depicted  in  Figure  18.  During  the  1st  step  the  last  node  (3  ),  calculates  its 
own  weight  w3  and  broadcasts  it  along  with  its  position  x3 .  In  the  2nd  step  the  2nd  node 

uses  the  received  information  to  solve  for  its  weight  w2 .  Then  it  broadcasts  >v2  and  its 
position  x2 .  The  procedure  goes  on  until  the  1st  node  calculates  its  own  weight  w1 
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Figure  18.  Last  phase  of  the  algorithm  to  implement  distributed  back  substitution  by 

the  sensor  nodes.  Bolded  and  underlined  xj ,  wt  above  the  nodes  denote  the 

positions  and  weights  that  are  already  stored  in  the  sensor.  Simple  xi,wi  denote 
the  broadcast  position  and  just  computed  weight. 


4.  Computational  and  Communication  Cost  Analysis 

The  computational  cost  of  the  procedure  presented  in  the  previous  section  can  be 
measured  in  number  of  instructions  needed  for  the  implementation.  It  is  known  that  the 
number  of  flops  for  the  solution  of  the  LS  problem  in  a  central  processor  is  given  by  [28] 

Nj  =2n2(m-n/3)  +  mn  +  n2  (71) 

where  the  first  term  is  for  the  determination  of  the  QR  factorization,  the  second  term  for 
the  update  of  the  desired  response  vector  with  the  Householder  transformations,  and  the 
last  term  accounts  for  the  back  substitution  to  solve  for  the  weights. 

The  significant  characteristic  of  the  proposed  distributed  algorithm  is  that,  as 
described  before,  no  redundant  computations  are  made.  The  Householder  transformations 

are  calculated  in  exactly  the  same  way  as  the  centralized  approach,  which  yields  the  same 
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number  of  computations.  In  the  centralized  approach,  however,  each  element  of  the 
distributed  stored  steering  matrix  D^n){6) ,  or  similarly  the  upper  triangular  matrix  Rx  of 

(67),  is  computed  by  a  single  sensor  and  no  redundant  computations  are  made.  For 
example,  in  the  distributed  approach,  the  first  node  does  not  construct  the  entire  matrix. 
In  the  first  phase,  it  uses  only  the  first  column  for  defining  the  first  Householder  matrix 
Hx ,  and  then  in  the  back  substitution  phase  it  calculates  only  the  first  row  of  the  Rx  using 

the  Hx .  Similarly,  the  second  node  calculates  H2  using  the  second  column  and  updates 
the  second  row,  which  is  used  for  the  back  substitution.  The  following  matrix  shows  the 
elements  of  Rx  with  respect  to  the  node  that  uses  and  calculates  them: 


1  1  1  •••  1 

0  2  2  •••  2 


0  0  0 


(72) 


Therefore,  no  additional  computations  are  made,  and  the  total  processing  cost  is  exactly 
equal  to  that  of  the  centralized  approach  given  in  (71).  Thus,  without  any  increase,  the 
processing  load  is  shared  among  the  nodes,  and  the  cluster  head  is  relieved  from  the 
heavy  computational  effort. 


The  number  of  instructions  Nt  can  be  used  to  express  the  processing  power  using 
the  definitions  of  Chapter  II  for  computational  cost.  Therefore,  from  (19)  and  (71),  the 
processing  power  Pp  is  given  by 


Pp  =  Nt  x  P.  =  ^2n2(m  —  n/3)  +  mn  +  n2~^  x  P.  .  (73) 


However,  this  distribution  of  the  processing  load  comes  with  a  tradeoff,  which  is 

an  increase  in  the  communication  load.  Note  that  throughout  this  work  the 

communication  cost  is  only  the  number  of  data  elements  that  need  to  be  transmitted  for 

the  implementation  of  an  algorithm.  The  transmissions  required  for  the  coordination  of 

the  sensors  in  the  distributed  algorithm  are  not  taken  into  account  for  the  calculation  of 

the  communication  cost.  In  the  centralized  approach,  the  cluster  head  collects  all  position 
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data  from  n  nodes  and  sends  back  n  weights,  so  the  number  of  transmitted  data 
elements  Nt  is  equal  to 

Nt=2n .  (74) 

On  the  other  hand,  in  the  distributed  algorithm,  the  number  of  transmitted  data 
elements  is  calculated  as  follows.  The  first  node  sends  the  H]  matrix,  or  equivalently  the 

mxl vector  vx  and  the  scalar  /?, ,  a  total  of  m  + 1  data  elements.  Similarly,  the  second 
node  broadcasts  the  matrix  H2  or  the  (m  - 1)  x  1  vector  y2  and  the  scalar  ,  a  total  of 
m  data  elements.  This  procedure  is  continued  until  the  (n  -  \)th  node,  which  transmits 
the  Hn  | .  In  general,  an  arbitrary  node  i  transmits  the  matrix  Hj ,  i.e.,  (m-i  + 1) 
elements  for  yf  and  one  for  /? .  The  total  number  of  transmitted  elements  by  all  nodes  for 
the  first  phase  of  the  QR  decomposition  is  given  by 


Cqr  =X(ffl-l  +  2)  =  (ffl  +  2)(n-l)-n^  ^  =  (m  +  2-n/2)(n-\)  .  (75) 

For  the  phase  of  the  back  substitution,  the  nth  node  transmits  its  position  and  its 
weight,  a  total  of  two  elements.  Similarly,  each  node  except  the  first  one,  broadcasts  its 
position  and  weight  to  the  other  nodes.  Thus,  the  total  number  transmitted  elements  for 
the  back  substitution  is  given  by 

CBs  =  2{n  ~  1)  •  (76) 

Finally,  the  total  number  of  data  elements  to  be  broadcast  during  the  implementation  of 
the  algorithm  is 

C  =  CqR  +  CBS  =(m  +  4-n/ 2 )(n  - 1)  .  (77) 

Using  again  the  definitions  from  Chapter  II  for  the  communication  power,  the 
power  consumption  due  to  communication  can  be  derived.  Assuming  that  each  element  is 
represented  by  b  bits,  the  number  of  transmitted  bits  is 
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(78) 


Ntb=Cxb 

and  from  (18)  the  communication  power  Pc  is  given  by 

Pc=NlbxPlb=(m  +  4-n/2)(n-l)xPtbxb.  (79) 

Therefore,  the  total  power  for  the  implementation  of  the  algorithm  is  the  sum  of  (73)  and 
(79) 

P  =  ^2n2  (m  -  n/3)  +  mn  +  n2~^x  P.  +  (m  +  4  -  n  /  2)(n  -  l)x  Ptb  xb  (80) 

and  depends  highly  on  the  specific  characteristics  Pt,Ptb  and  b  of  the  sensors.  Similarly, 
the  total  power  for  the  centralized  implementation  is  calculated  to  be 

P  =  ^2n2 (m  - n/3)  +  mn  +  «2  Jx P.  +  2nx Ptb  xb  .  (81) 

The  following  figures  summarize  the  results  and  compare  the  two  implementations. 

In  Figure  19,  the  number  of  required  instructions  AT  is  plotted  as  a  function  of  the 
number  of  approximation  points  m  for  a  10x1  array  of  sensors  while  in  Figure  20,  Nt  is 
plotted  as  a  function  of  the  number  of  sensors  n  ,  for  m  =  20  approximation  points.  The 
plots  for  the  processing  power  P  can  be  obtained  by  multiplying  the  number  of 

instructions  with  Pt  as  in  (73).  Another  important  point  is  that  the  main  processing  effort 

happens  during  the  QR  decomposition  phase  as  the  back  substitution  procedure  is  not  so 
demanding.  This  is  significant  because  the  QR  decomposition  is  computed  once  and  then 
the  results  are  kept  for  future  modifications.  For  example,  if  a  new  array  response  is 
required,  it  is  only  necessary  to  compute  the  updated  vector  cx  of  (65)  and  then  solve  the 
system  of  (66)  by  back  substitution;  Rx  does  not  change.  Therefore,  in  the  case  of 

communication  with  a  UAV,  which  is  moving  and  causing  the  desired  response  to  change 
continuously,  the  sensor  array  only  has  to  quickly  perform  the  back  substitution  phase  in 
order  to  compute  the  new  weights. 
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Figure  19.  Processing  cost  of  the  distributed  algorithm  as  a  function  of  the  number  of 

approximation  angles  and  for  n  =  1 0  sensors.  Multiplying  N i  by  P.  gives  the 
required  processing  power. 


Figure  20.  Processing  cost  of  the  distributed  algorithm  as  a  function  of  the  number  of 

sensors  and  for  fixed  number  of  approximation  angles  m  =  20 .  Multiplying  Ni 

by  P  gives  the  required  processing  power. 
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The  communication  cost  as  a  function  of  the  number  of  approximation  points  and 
the  number  of  sensors  is  depicted  in  Figures  21  and  22,  respectively.  Figure  21  shows  the 
communication  cost  for  a  10x1  sensor  array  in  terms  of  the  transmitted  elements. 
Assuming  that  each  data  element  is  represented  by  b  =  32  bits  in  a  single  precision 
floating  point  representation,  the  total  transmission  power  can  be  found  by  multiplying 
the  number  of  elements  C  by  Ptb*b  as  in  (78)  and  (79).  Similarly,  in  Figure  22,  the 

number  of  transmitted  elements  C  is  plotted  as  a  function  of  the  number  of  sensors  n 
for  m  =  20  .  Compared  to  the  centralized  approach,  the  communication  load  is  increased, 
especially  due  to  the  first  phase  of  the  QR  decomposition.  Flowever,  as  described  before 
for  the  processing  load,  this  phase  has  to  be  implemented  only  once.  Therefore,  any 
change  in  the  desired  response  needs  only  the  back  substitution  phase,  which  requires  a 
considerably  lower  number  of  elements  to  be  transmitted. 


Figure  21.  Communication  cost  of  the  distributed  algorithm  as  a  function  of  the 

number  of  approximation  angles  and  for  a  fixed  number  of  sensors  («  =  10 ). 
Multiplying  the  number  of  data  elements  with  Plbxb  gives  the  required 
transmission  power. 
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Figure  22.  Communication  cost  of  the  distributed  algorithm  as  a  function  of  the 

number  of  sensors  and  for  fixed  number  of  approximation  angles  ( m  =  20 ). 
Multiplying  the  number  of  data  elements  with  Ptbxb  gives  the  required 
transmission  power. 

In  Figures  23  and  24,  the  total  power  consumption  as  a  function  of  the  number  of 
approximation  points  and  the  number  of  sensors,  respectively,  is  plotted.  Assuming  that 
Pt  is  equal  to  one  unit  of  power,  and  using  the  ratio  Tjtp ,  the  factor  Ptb  can  be  substituted 

by  Ptb=7ltpy'I>i  ■  Dividing  the  total  power  P  of  (80)  by  P  yields  the  normalized  power 
Pn  in  terms  of  the  required  units  of  power  P  : 

Pn  =^  =  ^2n2(m-n/3)  +  mn  +  n2~^  +  (m  +  4-n/2)(n-\)xbxr]tp  (82) 

From  the  figures,  it  is  obvious  that  the  distributed  algorithm  requires  more  power 
than  the  centralized  one,  and  this  is  caused  by  the  additional  communication  load. 
However,  this  total  power  cost  is  shared  among  the  nodes  of  the  sensor  network  and  is 
not  undertaken  by  a  single  node.  For  example,  Figure  24  shows  that  a  deployment  of 
twenty  sensors  needs  about  five  times  more  power  for  the  distributed  approach  than  the 
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centralized  approach.  Therefore,  the  cluster  head  in  a  sensor  node  consumes  four  times 
more  power  in  the  centralized  scheme  than  the  average  sensor  nodes  consume  in  the 
distributed  one.  This  obviously  will  cause  the  cluster  head  to  fail  quickly,  which  means 
that  a  new  cluster  head  will  need  to  be  selected  and  all  the  computations  from  the 
beginning  will  need  to  be  re-done.  Thus,  the  increase  in  the  power  consumption  in  the 
distributed  approach  can  be  considered  reasonable  if  the  robustness  of  the  network  is  a 
crucial  requirement. 


Figure  23.  Normalized  power  Pn  (number  of  Pj )  for  both  distributed  and  centralized 

approaches  as  a  function  of  the  number  of  approximation  angles  for  a  fixed 
number  of  sensors  (n  =  20 ).  rj  is  assumed  to  be  200  and  b  =  32  bits. 
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Figure  24.  Normalized  power  Pn  (number  of  Pj )  for  both  distributed  and  centralized 

approaches  as  a  function  of  the  number  of  sensors  for  a  fixed  number  of 
approximation  angles  (m  =  20 ).  77  is  assumed  to  be  200  and  b  =  32 . 


B.  DISTRIBUTED  ITERATIVE  SCHEME  FOR  SOLVING  THE  LEAST 

SQUARE  PROBLEM 

In  this  selection,  a  distributed  scheme  based  on  an  iterative  solution  of  the  LS 
problem  is  presented  and  evaluated.  The  iterative  procedure  is  performed  in  steps  by  all 
the  nodes  in  the  sensor  array.  Starting  from  an  arbitrary  initialization  for  the  weights,  this 
method  quickly  converges  to  the  actual  solution  of  the  LS  problem. 

For  the  implementation  of  the  algorithm,  the  following  assumptions  are  made: 

a)  Each  node  can  calculate  its  position  from  a  reference  node  accurately. 

b)  Each  node  can  broadcast  information  to  other  nodes  in  the  cluster. 

c)  The  desired  array  response  for  a  set  of  directions  Qv...,0m  of  the  incoming 
signal  is  known. 
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d)  Errors  due  to  noise  during  the  communication  among  the  nodes  are  not  taken 
into  account. 

The  linear  system  of  (31),  which  is  to  be  solved  in  the  LS  sense,  is  considered  a 
minimization  problem  of  the  following  form: 

I  |2 

s  =  mm\Ax-b\  (83) 

x 

where 

A = D(0)H 

x  =  w  (84) 

b  =  Fd{6)  . 


1.  Proposed  Algorithm 

The  algorithm  is  based  on  various  parallel  methods  proposed  in  the  literature  for 
solving  the  LS  problem  [30],  [31],  [32],  [33], 

Let  A  be  an  in  x  n  matrix  while  x  and  b  are  n  x  1  and  m  x  1  vectors, 
respectively.  Each  column  of  A  is  denoted  4,  and  each  scalar  element  of  vector  x  is 
denoted  xt ,  for  /  =  I .  Then,  the  LS  problem  has  the  equivalent  form 

e  =  min^X;  +A2x2  +...  +  Anxn  -bf  .  (85) 


Suppose  that  xk)  is  an  approximation  to  the  solution  x  after  k  iterations  and  its 
elements  are  x\k) ,  for  i  =  Considering  an  arbitrary  element  x; ,  all  elements 

xlv..xM  are  updated  in  k  + 1  iterations  while  the  rest  xi,...xn  are  updated  in  k 
iterations.  Now,  (85)  can  be  written  as 


£  =  min 

v(*+l) 


l-l 


Edf 

7=1 


x<k+1)+A,x! 


.(*+i) 


+  Zd 

J=M 


7  7 


(86) 
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for  the  ( k  +  \)th  iteration  of  x\k+X) ,  which  gives  the  local  solution  for  xk+u .  The 
argument  of  the  term  on  the  right  hand  side  of  (86)  can  be  written  as 


/-I 


j= i 


j=i+ 1 


(87) 


-  Ak)) + + Z4/-v/Ai  -  * 

7=1 


J=l 


By  substituting  s^k)  as  the  step  or  correction 


„(*)  _  v(*+i)  _  v(*) 


(88) 


and  r(Av  1}  as  the  residual 


into  (86)  yields 


i-l 


=  Z4 

7=1 


j=i 


v(*) 

JXj 


8  =  mm 

„(*) 


XcW+r(i’M) 


(89) 


(90) 


Thus,  the  global  problem  of  finding  the  solution  for  x, , . . . ,  xn  in  (85)  is  equivalent 
to  solving  the  subproblems  of  (90),  which  can  be  assigned  to  the  corresponding  sensor 
nodes.  Indeed,  At  is  known  locally  to  the  sensor  node  i ,  sk)  is  the  locally  computed 
correction,  and  is  the  residual  after  the  i  —  \,h  node  has  completed  its  update  x^/ . 

The  residual  is  not  locally  available  in  the  node  i ;  however,  it  can  be  transmitted 

by  the  i-\th  sensor  node.  After  the  i th  node  calculates  the  new  solution  x\k\  it  sends 

the  updated  residual  rkJ)  to  the  other  nodes.  The  residual  can  be  shown  to  satisfy  the 
recursive  equation  [30] 


=  r 


(k,i- 1) 


( k ) 


(91) 


whereas  the  new  approximate  solution  is 


x,(*+1)=x,(i)+j,(t). 
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(92) 


Therefore,  the  i >h  node  is  assigned  a  column  .T,  of  the  matrix  A  and  assumes  an 
initial  solution  x.0) .  Assuming  that  r(0\  the  initial  estimation  for  the  residual  is  available, 
the  i th  node  solves  (90)  for  s\ 1 '  and  then  updates  its  solution  x\ 1 1  by  (92).  Following 
this  step,  the  updated  residual  r(U  is  sent  to  the  i  +  \th  node  in  order  to  update  its  solution 
x'+J .  The  procedure  continues  until  a  convergence  criterion  is  satisfied.  This  series  of 

approximations  converge  to  the  solution  x* ,  and  the  norm  of  the  residual  decreases 
continuously  [30],  [31].  This  procedure  is  summarized  in  Figure  25. 

a)  Divide  A  into  its  columns  At ,  for  i  =  l,...,n,  and  assign  one  to  each  node  i . 

b)  Initialize  x;<01 ,  for  i  =  \,...,n,  and  r(0) . 

c)  For  k  =  1  until  convergence 
For  i  =  1 

Solve  for  sik>  :  min  A  s(k>  +  . 

1  s(k)  ~ 

Update  the  residual :  r(k,l)  =  r(k,I~l)  +  Ats\k) . 

Send  the  updated  residual  to  all  nodes. 

Update  the  solution  xiu+l  )  =  x\k)  +  s(k) . 

Check  for  convergence. 

Figure  25.  Procedure  for  the  distributed  iterative  solution  of  the  LS  problem  (After 

Ref.  [30]). 

This  algorithm  is  highly  distributed,  and  it  can  be  implemented  in  a  sensor 
network  in  order  to  spread  the  computational  effort  equally  among  all  participating  sensor 
nodes. 


66 


From  the  notation  of  (84),  the  matrix  A  is  the  steering  matrix  D(6)H ,  each 
column  At  is  the  i"'  column  of  D(0)" ,  and  xt  is  the  weight  coefficient  wt.  From  (61), 

each  column  of  D(9)H  depends  only  on  the  position  of  the  node  and  the  set  of 
approximation  points.  This  column  is  available  locally  to  the  sensor  node.  Starting  from 
an  initial  estimation  for  the  residual  r(0) ,  the  first  node  can  solve  (90)  for  s,(1)  and  update 

its  own  weight  wf]  by  (92).  Then,  it  can  update  the  residual  r(U)  and  send  it  to  the  next 

node.  However,  since  the  residual  is  an  mx  1  vector,  each  transmission  of  the  residual 
needs  the  transmission  of  m  elements,  which  is  a  considerable  amount  of  communication 
load.  Another  approach  for  the  transmission  of  the  residual  seems  more  efficient.  Assume 
that  all  sensors’  positions  are  initially  broadcast  to  all  sensor  nodes  so  that  each  node 
knows  the  position  of  each  other  node;  in  this  way,  they  can  construct  any  of  the  columns 
of  the  steering  matrix.  Then  the  ith  sensor  node  does  not  need  to  send  the  entire  residual 
rkj)  but  only  the  scalar  correction  s\k] ,  and  all  other  nodes  can  reconstruct  the  residual 
by  repeating  the  (91).  For  example,  the  second  node  is  not  required  to  send  the  residual 
rk2)  but  just  the  scalar  correction  s(2k} ,  and  then  the  third  node  will  reconstruct  the 
residual  by  computing  (91): 


Jk,  2) 


=  r(k,l)  +  A2s2 


(k) 


(93) 


Therefore,  in  each  iteration  step,  only  one  scalar  element  has  to  be  broadcast,  and 
the  rest  of  the  computations  are  performed  locally.  This  distributed  determination  of  the 
weights  in  a  sensor  network  is  summarized  in  Figure  26. 
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a)  All  nodes  broadcast  their  position;  each  node  can  construct  all  columns 
D, ,  for  i  =  l,...,n  ,  of  D(6)H  . 

b)  Initialize  w f0) ,  i  =  \,...,n  and  r(0) . 

c)  For  k  =  1  until  convergence 

For  i  =  1 

Receive  s(kJ  from  i  - 1  node. 

Reconstruct  the  residual:  r(k,‘~l)  =  rkJ  21  +  D,  . 

Solve  for  s' k '  :  min  Ds\ k '  +  r'kJ  "  . 

Update  the  residual:  r(k,l)  =  rkJ  "  +  Dt :s\k) . 

Send  the  only  the  correction  sik)  to  all  nodes. 

Update  the  solution  w{k+']  =  w'k)  +  s(k  ) . 

Check  for  convergence. 

Figure  26.  Procedure  for  the  proposed  distributed  iterative  solution  of  the  LS  problem 

in  a  WSN  environment. 

2.  Computational  and  Communication  Costs 

The  performance  of  the  distributed  iterative  procedure  for  the  computation  of  the 
weights  in  a  random  sensor  array  is  shown  in  Figures  27  and  28.  In  Figure  27,  the 
residual  norm  is  plotted  for  each  local  iteration  for  a  10x1  array  of  sensors.  It  is  obvious 
that  the  residual  norm  converges  to  the  actual  residual  after  about  thirty  local  iterations  or 
about  three  complete  iterations;  one  complete  iteration  is  finished  when  all  ten  nodes 
perform  a  local  iteration.  Similar  results  are  plotted  in  Figure  28  where  the  norm  of  the 
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error  e(k )  between  the  approximate  w  k>  and  the  actual  solution  w  calculated  after 
each  complete  iteration  k  is 


e{k)  = 


(94) 


Figure  27.  Convergence  of  the  residual  norm  to  the  actual  residual,  indicating  that  the 

algorithm  converges  to  the  real  solution.  After  3  complete  iterations  (30  local)  the 
residual  has  converged. 
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Figure  28.  Convergence  of  the  norm  of  the  error  e( k)  between  the  approximate  wk> 

and  actual  solution  yv*  calculated  after  each  complete  iteration  k ,  indicating  that 
the  algorithm  converges  to  the  real  solution. 

The  computational  cost  can  be  measured  in  terms  of  the  number  of  instructions 
needed  for  the  implementation.  Assuming  that  each  node  solves  its  local  LS  problem  with 
a  QR  factorization  obtained  by  Householder  transformations,  the  cost  per  sensor  N  is 

given  by  (71),  which  for  n  =  1  yields 

Nps  =  2{m  - 1/3)  +  m  +  \  .  (95) 

Note  that  the  last  two  terms  stand  for  the  backsubstitution,  so  they  are  performed  in  each 
iteration,  thus  (95)  becomes 

Nps=2(m-l/3)  +  k(m  +  l).  (96) 

Furthermore,  each  node  has  to  reconstruct  the  in  x  I  residual  vector,  so  an  additional 
number  of  2m  instructions  are  needed.  Finally,  the  number  of  instructions  per  sensor 
results  in 


70 


Nps  =  2(m  - 1/3)  +  k(3m  + 1)  . 


(97) 


Therefore,  the  total  processing  cost  for  all  the  nodes  Nt  is  equal  to 

Ni=nx  Nps  =  n  [2(m  - 1/3)  +  k(3  m  + 1)]  (98) 

while  the  total  processing  power  is 

Pp=NjXPi=n  [2 (m  - 1/3)  +  k(3m  + 1)]  x  .  (99) 

Using  the  above  equation,  the  processing  cost  in  number  of  instructions  for  the 
implementation  of  the  algorithm  is  plotted  as  a  function  of  the  number  of  sensors  n ,  the 
number  of  approximation  angles  m,  and  the  number  of  iterations  k .  In  Figure  29,  the 
number  of  instructions  for  the  distributed  approach  become  fewer  than  the  centralized 
approach  after  a  certain  point,  which  is  m  =  15;  the  simulation  scenario  is  for  a  10x1 
array  of  sensors  and  the  cost  has  been  computed  for  k  =  5  iterations.  Similarly,  in  Figure 
30,  the  processing  cost  for  the  centralized  approach  grows  larger  than  the  distributed 
approach  as  the  number  of  the  sensors  increase  (m  =  20  and  k  =  5). 
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Figure  29.  Processing  cost  of  the  distributed  algorithm  as  a  function  of  the  number  of 

approximation  angles  for  n  =  10  sensors  and  k  =  5  iterations.  Multiplying  the 
number  of  instructions  by  P.  gives  the  required  processing  power. 


Figure  30.  Processing  cost  of  the  distributed  algorithm  as  a  function  of  the  number  of 

sensors  for  m  =  20  approximation  angles  and  k  =  5  iterations.  Multiplying  the 
number  of  instructions  by  Pt  gives  the  required  processing  power. 

72 


Figure  3 1 .  Processing  cost  of  the  distributed  algorithm  as  a  function  of  the  number  of 

iterations  for  n  =  10  sensors  and  m  =  20  approximation  angles.  Multiplying  the 
number  of  instructions  by  P  gives  the  required  processing  power. 

Finally,  Figure  31  shows  the  increase  of  the  processing  cost  with  the  number  of 
iterations,  which  indicates  that  the  cost  for  the  distributed  approach  grows  higher  than  the 
cost  for  the  centralized  one  after  a  specific  number  of  iterations  (&  =  5in  this  case);  a 
10x1  sensor  array  is  deployed  and  twenty  approximation  points  are  used. 

The  communication  cost  for  the  distributed  approach  can  be  derived  as  follows. 
Initially,  each  node  broadcasts  its  own  position,  so  a  total  of  n  data  elements  are 
transmitted.  Then,  after  a  local  iteration  step  is  finished,  the  scalar  correction  sjk)  has  to 

be  sent,  so  for  k  complete  iterations,  kn  elements  are  transmitted.  Therefore,  the  number 
of  transmitted  elements  is 

C  =  (k  +  l)n  (100) 

and  it  does  not  depend  on  the  number  of  approximation  angles  m  .  Assuming  that  each 
element  is  represented  by  b  bits,  the  total  transmission  power  Pc  is  given  by 
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Pc  =  C  x  Ptb  x  B  =  {k  +  \)n  x  Ptb  x  b . 


(101) 


In  Figure  32,  the  communication  cost  is  plotted  as  a  function  of  the  number  of 
sensors.  The  increased  communication  load  for  the  distributed  approach  is  a  reasonable 
tradeoff  considering  that  this  distributed  approach  relieves  the  central  processing  node 
from  the  entire  computational  load.  Similarly,  Figure  33  shows  the  effect  of  the  number 
of  iterations  on  the  communication  cost.  Obviously,  the  speed  of  convergence  positively 
affects  the  reduction  of  the  transmitted  elements. 

The  total  power  for  the  implementation  of  the  algorithm  is  the  sum  of  the 
processing  power  given  by  (99)  and  the  transmission  power  given  by  (101): 

P  =  n  [2  (m  - 1/3)  +  k(3m  +  l)]xP.+(k  +  l)nxPtbxb  (102) 

Considering  the  ratio  77  and  substituting  Ptb  from  (21),  the  normalized  power  Pn 
as  a  number  of  required  units  of  power  P.  is  given  by 

Pn  =  =  n[2(m -1/3)  +  k(3m  +  l)]xP.  +(k  +  l)nxrjtp  xb  .  (103) 
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Figure  32.  Communication  Cost  of  the  distributed  algorithm  as  a  function  of  the 

number  of  sensors  for  m  =  20  approximation  angles  and  k  =  5  iterations. 
Multiplying  the  number  of  data  elements  with  Ptb  x  b  gives  the  required 
transmission  power. 


Number  of  iterations 


Figure  33.  Communication  cost  of  the  distributed  algorithm  as  a  function  of  the 

number  of  iterations  for  « =  1 0  sensors  and  m  =  20  approximation  angles. 
Multiplying  the  number  of  data  elements  with  Plbxb  gives  the  required 
transmission  power. 
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In  Figures  34-36,  the  total  normalized  power  needed  for  the  implementation  of  the 
algorithm  is  plotted  as  a  function  of  the  number  of  approximation  angles,  the  number  of 
sensors  and  the  number  of  iterations,  respectively.  In  all  cases,  the  power  consumption  of 
the  distributed  algorithm  is  larger  than  that  of  the  centralized  approach.  However,  the 
important  characteristic,  as  in  the  previous  algorithm,  is  that  the  total  amount  of  power  is 
shared  among  the  nodes;  consequently,  the  cluster  head  is  relieved  of  the  computational 
burden.  For  instance,  in  Figure  35,  where  m  =  20  and  k  =  5  ,  the  power  for  the  distributed 
approach  is  almost  three  times  larger  than  the  centralized  approach,  but  this  is  shared  by 
n  =  20  nodes.  Thus,  the  average  power  consumption  of  an  arbitrary  node  in  the 
distributed  network  is  six  times  less  compared  to  the  cluster  head’s  consumption  in  the 
centralized  architecture.  In  a  centralized  approach,  there  is  a  considerably  higher 
possibility  that  the  cluster  head  will  exhaust  its  limited  battery  and  then  the  beamforming 
problem  will  have  to  be  computed  from  scratch.  In  summary,  the  increased  power 
consumption  is  reasonable  enough  to  consider  the  distributed  approach  a  viable  solution 
if  robustness  is  required  in  the  network. 
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Figure  34.  Normalized  power  Pn  (number  of/>)  for  both  distributed  and  centralized 
approaches  as  a  function  of  the  number  of  approximation  angles  for 
n  =  10  sensors  and  k  =  5  iterations.  rjtp  is  assumed  to  be  200  and  b  =  32  bits. 
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Figure  35.  Normalized  power  Pn  (number  of  Pt)  for  both  distributed  and  centralized 
approaches  as  a  function  of  the  number  of  sensors  for  m  =  20  approximation 
points  and  k  =  5  iterations,  rj  is  assumed  to  be  200  and  b  =  32  bits. 


Figure  36.  Normalized  power  Pn  (number  of  P.)  for  both  distributed  and  centralized 
approaches  as  a  function  of  the  number  iterations  for  n  =  10  sensors  and 
m=  20  approximation  points,  rj  is  assumed  to  be  200  and  b  =  32  bits. 
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C.  SUMMARY 

In  this  chapter,  two  distributed  approaches  were  proposed  for  the  solution  of  the 
LS  problem  of  the  array  weight  computation.  The  presented  algorithms  were  evaluated  in 
terms  of  the  processing  cost,  the  communication  cost  and  the  total  power  consumption 
and  then  compared  to  the  centralized  implementation. 

In  the  first  approach  (distributed  QR  decomposition),  the  processing  cost  is  the 
same  as  in  the  centralized  approach,  but  there  is  a  tradeoff  of  increased  communication 
effort.  However,  only  the  first  phase  of  the  algorithm  is  power  demanding,  and  potential 
modifications  to  the  desired  response  do  not  require  much  power  for  the  recalculation  of 
the  weights. 

In  the  second  approach  (iterative  solution),  there  is  rapid  convergence  to  the 
actual  solution,  which  yields  a  reduction  of  the  processing  cost  if  the  number  of  iterations 
is  kept  low.  The  simulation  results  show  that  only  3  or  4  iterations  are  needed  for  the 
convergence  of  the  algorithm,  which  results  in  considerably  lower  processing  power. 
However,  the  communication  cost  is  still  higher  when  compared  to  the  centralized 
approach. 

To  sum  it  up,  these  two  approaches  require  significantly  lower  average  power  per 
sensor  node  by  efficiently  sharing  the  power  consumption  among  the  nodes. 
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V.  CONCLUSIONS 


The  operational  scenario  adopted  in  this  work  assumes  that  a  number  of  sensor 
nodes  are  randomly  deployed  in  an  area  of  interest  in  order  to  collect  information  about 
various  kinds  of  objects.  The  acquired  data  has  to  be  collected  by  an  oveflying  UAV; 
however,  single  sensors  do  not  have  sufficient  power  capabilities  in  order  to  establish 
communication  with  the  UAV.  Therefore,  they  are  organized  into  clusters  and  cooperate 
in  order  to  function  as  an  array  of  sensor  nodes. 

The  effect  of  position  errors  on  the  performance  of  the  random  sensor  array  was 
analyzed,  and  the  need  for  beamforming  techniques  which  effectively  mitigate  that  effect 
was  discussed.  Since  reliability  and  robustness  in  a  sensor  network  environment  are 
crucial,  two  distributed  algorithms  for  beamforming  that  efficiently  manage  to  share  the 
processing  load  among  the  sensor  nodes,  compared  to  centralized  approaches  which 
assign  the  entire  effort  to  a  single  node,  were  presented.  A  simulation  model  was  created 
and  implemented  in  the  MATLAB  environment  to  evaluate  the  performance  of  the 
proposed  algorithms. 


A.  SIGNIFICANT  RESULTS 

The  simulations  showed  that  the  sidelobes  in  the  array  response  increase  as  a 
function  of  the  “randomness”  of  the  sensor  array.  Thus,  as  the  mean  deviation  from  the 
uniform  array  was  increased,  the  average  sidelobe  magnitudes  also  increased.  These 
results  validate  the  theoretical  results  of  random  arrays  found  in  the  literature.  Another 
important  point  is  that  for  the  LS  beamformer,  a  subset  of  approximation  points  can  yield 
almost  the  same  solution  as  larger  sets.  This  offers  significant  reduction  to  the  required 
processing  and  transmission  power,  which  are  crucial  in  sensor  networks. 

Based  on  the  performance  analysis,  the  two  proposed  distributed  algorithms  can 

effectively  share  the  processing  load  among  the  nodes.  The  first,  a  distributed 

implementation  of  the  QR  decomposition,  has  the  same  processing  cost  as  the  centralized 

one.  The  second  approach,  based  on  an  iterative  method  of  computing  the  weight  vector 
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in  the  LS  sense,  converges  quickly  to  the  actual  solution  and  achieves  reduction  of  the 
total  processing  cost  compared  to  that  of  the  centralized  one. 

For  both  algorithms,  the  tradeoff  is  the  increased  transmission  power,  causing  an 
overall  increase  in  the  total  power  consumption  in  the  network.  This  total  power, 
however,  is  shared  among  the  sensor  nodes;  therefore,  the  average  power  needed  by  a 
sensor  node  in  the  distributed  implementation  is  lower  than  the  power  needed  by  the 
cluster  head  in  the  centralized  approach.  Consequently,  the  network’s  susceptibility  to 
failures  due  to  excessive  power  consumption  is  greatly  reduced. 


B.  FUTURE  WORK 

Throughout  this  work,  several  assumptions  were  made,  such  as  the  nodes  can 
compute  their  positions  without  errors,  and  the  communication  between  them  is  not 
affected  by  noise.  A  future  effort  may  examine  the  effect  of  these  errors  on  the 
computation  of  the  weight  vector  and  consequently  on  the  array  performance. 

In  this  work,  the  set  of  approximation  points  were  selected  based  on  uniform 
sampling,  but  there  are  other  choices,  such  as  using  a  non-uniform  grid,  which  may 
require  fewer  points  with  similar  performance.  Initial  results  showed  that  certain 
approximation  points,  which  have  the  physical  meaning  of  direction  angles,  may  be  more 
important  for  the  approximation  of  the  desired  response  than  others.  These  issues  may  be 
further  investigated  in  a  future  study. 

The  topology  of  WSN  changes  dynamically  due  to  frequent  additions  and 
withdrawals  of  sensor  nodes;  some  of  them  may  switch  in  or  off  sleep  mode  and  some 
other  may  fail  because  of  the  harsh  environmental  conditions  or  because  of  exhausted 
battery.  For  these  scenarios,  the  processing  and  communication  cost  for  the  update  of  the 
weight  vectors  can  be  investigated.  Also  for  the  distributed  algorithms,  there  is  no  need  to 
solve  the  beamforming  problem  from  scratch;  if  the  array  topology  changes,  it  is 
important  to  analyze  the  effects  on  the  costs  and  consequently  the  power  needed  for 
modifying  the  weight  vector  after  a  node  has  added  or  withdrawn  from  the  array. 
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In  this  work,  the  emphasis  was  given  to  the  data  independent  beamforming 
techniques,  such  as  the  LS  approximation  of  the  desired  response.  There  are  many  proven 
data  dependent  techniques  for  which  the  weight  vector  can  be  determined  adaptively 
[11],  [12]  and  a  distributed  implementation  of  these  methods  could  be  examined. 

An  array  of  M  elements  can  be  used  to  form  a  beampattem  with  exactly  M  - 1 
narrowband  nulls.  These  nulls  can  be  placed  towards  the  directions  of  incoming 
interferences  in  order  to  suppress  them.  Since  there  are  straightforward  polynomial  based 
techniques  for  placing  those  nulls  where  desired  in  the  array  space,  it  would  be 
interesting  to  investigate  them  in  a  future  work. 

Finally,  the  communication  cost  was  defined  as  a  function  of  the  data  elements 
that  need  to  be  transmitted  for  the  implementation  of  the  algorithms.  However,  there  is 
also  a  networking  cost  which  consists  of  parameters  such  as  the  packet  overhead  and  the 
retransmissions  due  to  collisions  and  errors,  which  all  add  to  the  power  consumption. 
Since  the  implementation  of  a  distributed  beamforming  algorithm  may  be  prohibitive  by 
a  high  networking  cost,  it  will  be  interesting  to  investigate  the  effect  of  the  networking 
cost  to  the  overall  power  consumption. 
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APPENDIX 


MATLAB  SOURCE  CODE 


This  appendix  lists  all  MATLAB  programs  used  in  this  work 

•  Array2D.m : 


Array2D . m 

Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  file  generates  the  array  beampattern  for  an 
array  with  randomly  positioned  elements 
The  weights  are  computed  using  the  LS  approach  of 
the  desired  response 


clear  all 
close  all 
clc 


Filename : 
Author : 


Date : 

Description : 


PARAMETERS 


oooooooooooooooooooooooooooooooooooooooooooooooo 


global  c  f  1  b  Im  Nx  Ny 
c=3e8 ; 


f=2e9 ; 
l=c/f; 
b=2*pi/l; 


GdBavg=zeros (181,181) ; 
GdBerravg=zeros  (181, 181)  ; 
GdBerrlinavg=zeros (181, 181) ; 
GdBref avg=zeros (181, 181) ; 
GdBlsavg= zeros (181,181)  ; 
GdBLSavg=zeros (181, 181)  ; 
GdBiteravg=zeros (181, 181, 16)  ; 


GdBLSlavg=zeros (181, 181)  ; 
GdBLS2avg=zeros (181, 181)  ; 
GdBLS3avg=zeros (181, 181)  ; 
GdBLS4avg=zeros (181, 181)  ; 
GdBLS5avg=zeros (181, 181)  ; 


wwls= [ ] ; 
wwnu= [ ] ; 
WW1=  [  ]  ; 
WW2= [ ] ; 
WW3= [ ] ; 
WW4= [ ] ; 
WW5=  [  ]  ; 
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o 

o 


9'2'9'2'9'2'9'2'9'2'9' 

ooooooooooo 


END  OF  PARAMETERS 


ooooooooooo 


o 

o 


oooooooooooooooooooooooooooooooooooooooooooo 
9'9'9'9'9'9'9'9'S'  TMDTTT1  p*  u  n  T  P1  T71  Q 

ooooooooo  -L  IN  Jr  U  ±  ^nW±L^JliO  ooooooooooooooo 


o 

o 


P=inputparameters ; 


Nx=P  (1) ; 
Ny=P  (2)  ; 


Nx  number 
Ny  number 


of  array 
of  array 


elements 

elements 


in  x  direction 
in  y  direction 


XG=P ( 3 ) 

r 

o 

o 

linear. 

(2) :  From 

xe=P ( 4 ) 

; 

o 

o 

perfect 

linear) 

in 

xe=xe 

i— 1 

o 

o 

ye=P ( 5 ) 

r 

o 

o 

perfect 

linear) 

in 

ye=ye/100 ; 


Generation  of  position 
cratch 

Percentage  of  position 
x  direction  (%) 

Percentage  of  position 
y  direction  (%) 


(1) : Deviation  from  perfect 
error  (with  respect  to 

error  (with  respect  to 


xest=P(6);  %  Percentage  of  estimated  position  error 

to  actual)  in  x  direction  (%) 
xest=xest/100 ; 

yest=P(7);  %  Percentage  of  estimated  position  error 

to  actual)  in  y  direction  (%) 
yest=yest/100 ; 


(with  respect 

(with  respect 


thetaO=P(8);  %  Elevation  angle  theta  (degrees) 

theta0=theta0*pi/180 ; 

phiO=P(9);  %  Azimuth  angle  phi  (degrees) 

phi0=phi0*pi/180; 

phi_ang=P ( 10 ) ;  %  Angle  phi  for  beampattern 


te=P (11) ; 
pe=P ( 12 ) ; 


%  Angle  error  in  theta  (+-  degrees) 
% ('Angle  error  in  phi  (+-degrees) 


NumIter=P (13) ; 


%Number  of  iterations  (for  average  beampattern) 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 


ooooooooooooo 


END  OF  INPUT  CHOICES 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 


oooooooooooo 


ooooooooooooooooooooooooooooo 


o  o  o  o 
o  o  o  o 


o 

o 


o 

o 


o 

o 


o 

o 


o 

o 


o 

o 


oooooooo 


POSITION  GENERATION 


ooooooooooooooooooooooooooooooooooooooooooooooooooo 
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%%%%%%%%%  Uniform  array  (reference)  %%%%%%%%%% 

dx=l/2;  %  ideal  distance  lamda/2  in  x-direction 
xn= ( 0 : Nx-1 ) *dx; 
xn=repmat (xn',1, Ny) ; 
xn=reshape (xn, Nx*Ny, 1) ; 

dy=l/2;  %  ideal  distance  lamda/2  in  y-direction 

yn= ( 0 : Ny-1 ) *dy  ; 
yn=r epma t ( yn , Nx  ,  1 )  ; 
yn=reshape (yn,Nx*Ny,  1)  ; 


ooooooooo 


End  of  uniform  array 


ooooooooooo 


for  NI=1 : Numlter ; 

%%%%%%%%  DISTANCE  DEVIATION  %%%%%%%%% 
if  XG==1 ; 

devx=xe*l* ( rand (Nx*Ny, 1 ) -0 . 5 ) ;  %  Random  deviation 

-lamda/2  to  xe%  lamda/2 

x=xn+devx;  %  Real  positions  in  x-direction 

devy=ye*l* (rand (Nx*Ny, 1) -0 . 5) ;  %  Random  deviation 

lamda/2  to  lamda/2 

y=yn+devy;  %  Real  positions  in  x-direction 

elseif  XG==2;%  2nd  option  -  Firstly  construct  x,y 
assume  linear 

[x, y] =rand_inter_dist (Nx, Ny) ; 

end 

%%%%  Move  reference  node  to  the  axes  center  %%%%%%% 

x=x-x ( 1 ) ;  % 

y=y-y(l);  % 


%%%%%%  Estimated  position  with  defined  error  %%%%%%% 

xerr=x . * ( 1+xest* ( 2 *rand (Nx^Ny, 1 ) -1 ) ) ;  %  error  position  in  x 

direction  with  respect  to  actual 

yerr=y . * ( 1+yest* (2^rand (Nx*Ny, 1 ) -1 ) ) ;  %  error  position  in  y- 

direction  with  respect  to  actual 


End  of 


Estimated  position  with  defined 
END  OF  DISTANCE  DEVIATION 


error 


o  o  o  o  o  o 
o  o  o  o  o  o 


ooooooooooo 


from  xe% 


from  - 


and  then 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%  WEIGHTS  %%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%  Amplitudes 


0.0.0. 
o  o  o 


Im=ones (Nx*Ny, 1) ; 


%  Amplitudes 


wm=weights2 (x, y, thetaO , phiO ) ;  %  correct  weights 

wmerr=weights2 (xerr ,  yerr ,  thetaO , phiO ) ;  %  wrong  weights,  deviation 

from  actual  position  :  xe% 

wmerrlin=weights2 (xn, yn, thetaO , phiO ) ;  %  wrong  weights,  assume 

perfect  linear 

wref =weights2 (xn, yn, thetaO , phiO ) ;  %  Reference  weights  (uniform  array) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%  END  OF  WEIGHTS  %%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%  GAIN  %%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


G=gain2D (wm, x, y) ;  %  Gain  from  correct  weights 

Gerr=gain2D (wmerr , x, y) ;  %  Gain  from  wrong  weights,  deviation  from 

actual  position 

Gerrlin=gain2D (wmerrlin, x, y) ;  %  Gain  from  wrong  weights,  assumed 

perfect  linear 
Gref=gain2D (wref, xn, yn) ; 


S'9'9'9'9'9'9'9'  n  -a  \  ra  9-  9-  o,  o.  o,  o.  o,  o.  o,  o.  o. 

oooooooo  balll  ^  UU_)  )  ooooooooooo 

GdB=10*logl0 (G/max (max (G) ) ) ;  %  Gain  from  correct  weights  (dB) 

GdBerr=10^1ogl0 (Gerr/max (max (G) ) ) ;  %  Gain  from  wrong  weigths, 

deviation  from  actual  position  (dB) 

GdBerrlin=10*logl0 (Gerrlin/max (max (G) ) ) ;  %  Gain  from  wrong  weights, 
assumed  perfect  linear  (dB) 

GdBref=10*logl0 (Gref /max (max (G) )) ;  %  Gain  for  reference  (uniform 

array) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%  END  OF  GAIN  %%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%%%%%%%%  LS  ESTIMATION  OF  THE  WEIGHTS,  GIVEN  DESIRED  RESPONSE 

ooooooooooo 


theta=-90 : 90 ; 
th=theta*pi/180  ; 


dn=exp ( j  *b* (xn*sin ( th) *cos (phiO ) +yn*sin ( th)  *sin  (phiO )  )  )  ; 
steering  vector  for  ULA 


Fdes=wref 1 *dn; 


d=exp (j*b* (x*sin (th) *cos (phiO ) +y*sin ( th) *sin (phiO ) ) ) ;  % 

steering  vector 

ww=Fdes/ d; 

ww=ww';  %  or  ww=inv (d*d ' ) *d*Fdes ’ ; 


wwls= [wwls  ww]  ; 

Gls=gain2D (ww, x, y) ; 
GdBls=10*logl0 (Gls/max (max (G) ) ) ; 


thO=P (8) ; 

%%%%%%%  LS  with  fewer  approximation  points  %%%% 

%%  Uniform  spacing  %% 

0,0  /  -i  \  oo 

o  o  \  ±  )  o  o 

dt=2  ; 

rla=thO : -dt : —  9  0 ; 
rla=f lipdim (rla, 2 ) ; 
rlb=thO+dt : dt : 90 ; 
rl= [rla  rib] ; 

tl«rl*pi/180; 

DNl=exp ( j  *b* (xn*sin ( tl ) *cos (phiO ) +yn*sin ( tl ) *sin (phiO ) ) ) 
FDESl=wref 1 *DN1 ; 

Dl=exp (j*b* (x*sin (tl) *cos (phiO ) +y*sin ( tl ) *sin (phiO ) ) ) ; 

wwl=FDESl/Dl; 
wwl=wwl 1 ; 

WW1=[WW1  wwl] ; 
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GLSl=gain2D (wwl , x, y) ; 
GdBLSl=10*logl0 (GLSl/max (max (G) ) )  ; 


0,0  /  O  \  0,0 

o  o  \  Z  ;  o  o 

dt=4  ; 

r2a=th0 : -dt : —  9  0 ; 
r2a=f lipdim ( r2a, 2 ) ; 
r2b=th0+dt : dt : 90/ 
r2= [ r2a  r2b] ; 

t2=r2*pi/180 ; 

DN2=exp ( j  *b* (xn*sin ( 1 2 ) *cos (phiO ) +yn*sin ( 1 2 ) *sin (phiO ) ) ) 
FDES2=wref 1  * DN 2 ; 


D2=exp ( j  *b* (x*sin ( 1 2 ) *cos (phiO ) +y*sin ( t2 ) *sin (phiO ) ) ) ; 

ww2=FDES2/D2; 
ww2=ww2 1 ; 

WW2=[WW2  ww2 ] ; 

GLS2=gain2D (ww2  ,  x,  y)  ; 

GdBLS2=10*logl0 (GLS2/max (max (G) ) ) ; 

0,0,  /  q  \  oo, 

o  o  ^  ;  o  o 

dt=6 ; 

r3a=th0 : -dt : —  9  0 ; 
r3a=f lipdim ( r3a, 2 ) ; 
r3b=th0+dt : dt : 90/ 
r3= [ r3a  r3b] ; 

t3=r3*pi/180 ; 

DN3=exp (j*b* (xn*sin ( 1 3 ) *cos (phiO ) +yn*sin ( 1 3 ) *sin (phiO ) ) ) 
FDES3=wref 1  * DN  3 ; 


D3=exp ( j  *b* (x*sin (t3) *cos (phiO) +y*sin (t3) *sin (phiO) ) ) ; 

ww3=FDES3/D3; 
ww3=ww3 1 ; 

WW3=[WW3  ww3 ] ; 

GLS3=gain2D (ww3 , x, y) ; 

GdBLS3=10*logl0 (GLS3/max (max (G) ) ) ; 
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dt=8 ; 

r4a=th0 : -dt : —  9  0 ; 
r4a=f lipdim ( r4a, 2 ) ; 
r4b=th0+dt : dt : 90/ 
r4= [ r4a  r4b] ; 

t4=r4*pi/180 ; 

DN4=exp ( j  *b* (xn*sin ( 1 4 ) *cos (phiO ) +yn* sin ( 1 4 ) *sin (phiO ) ) ) 
FDES4=wref 1 *DN4 ; 

D4=exp (j*b* (x*sin ( 1 4 ) *cos (phiO ) +y*sin ( 1 4 ) *sin (phiO ) ) ) ; 

ww4=FDES4/D4; 
ww4=ww4 1 ; 

WW4  = [WW4  ww4 ] ; 

GLS4=gain2D (ww4  ,  x,  y)  ; 

GdBLS4=10*logl0 (GLS4/max (max (G) ) ) ; 


9-  9-  (  R  \  0,0, 

o  o  [  )  o  o 

dt=10 ; 

r5a=th0 : -dt : —  9  0 ; 
r5a=f lipdim ( r5a, 2 ) ; 
r5b=th0+dt : dt : 90/ 
r5= [ r5a  r5b] ; 

t5=r5*pi/ 180 ; 

DN5=exp ( j  *b* (xn*sin ( 1 5 ) *cos (phiO) +yn*sin ( 1 5 ) *sin (phiO) ) ) 
FDES5=wref 1 *DN5 ; 

D5=exp ( j  *b* (x*sin ( 1 5 ) *cos (phiO) +y*sin (t5) *sin (phiO) ) ) ; 

ww5=FDES5/D5 ; 
ww5=ww5 ' ; 

WW5  = [WW5  ww5 ] ; 

GLS5=gain2D (ww5 , x, y) ; 

GdBLS5=10*logl0 (GLS5/max (max (G) ) ) ; 


%%%  Non  uniform  %%%% 
thetal=th0-12 : 4:th0+12; 
theta2=-90 : 15 : 90 ; 
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f t= (theta2<th0-12 ) | ( theta2>th0+12 ) ; 
thetaf= [ thetal  theta2  (ft) ] ; 

THETA=sort (thetaf )  ; 

TH=THETA*pi/ 180 ; 

DN=exp ( j  *b* (xn*sin (TH) *cos (phiO ) +yn*sin (TH) *sin (phiO ) ) ) 
FDES=wref 1 *DN; 

D=exp ( j  *b* (x*sin (TH) *cos (phiO ) +y*sin (TH) *sin (phiO ) ) ) ; 

WW=FDES/D; 

WW=WW ' ; 


wwnu=[wwnu  WW]  ; 


GLS=gain2D (WW, x, y) ; 
GdBLS=10*logl0 (GLS/max (max (G) ) ) ; 


ooooooooooooooooooo 


END  OF  LS 


oooooooooooooooooooooo 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


S'9'9'9'9'9'9'9'9'9'9'9'9'  PIOTC  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9- 

ooooooooooooo  l:  LU  ±  Q  ooooooooooooo 


figure  ( 10 )  ; 

plot (xn, yn, ’o') 
hold  on; 
plot (x, y, '  rx '  ) 
hold  on; 

plot (xerr , yerr , 'mp') 
grid  on 
axis  equal 

title('Fig.l  Sensor  array','Fontsize',12); 

legend (' Perfect  linear Actual  Position Wrongly  estimated’); 


figure (20) ; 

plot (theta, GdB ( : , phi_ang+90+l ) , ’ Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBerr ( : , phi_ang+90+l ) , ' r- . ’ , ’ Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBerrlin ( : , phi_ang+90+l ) , ' m : ’ , ' Linewidth ’ , 2 ) ; 
hold  on; 

plot (theta, GdBref ( : , phi_ang+90+l ) , ’ g- ' , ’ Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBls ( : , phi_ang+90  +  l ) ,  ' c- ’  ,  ' Linewidth ' , 2 ) ; 
axis ( [-85  85  -50  5] ) ; 
grid  on 

title ('Fig. 2  :  Beampattern  for  N  linear  array  elements  and  given 
\phi ' , ' Fonts ize ' , 12 ) ; 
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xlabel ( ' \ theta  (degrees )  ’  ,  ’ Font size '  ,  12 ) ; 
ylabel (' Power  Gain  (dB) ,  '  ,' Fontsize 1 , 12 ) ; 

legend (' Correct ' ,  ' Wrongly  estimated Assumed  perfect  linear ' ,  ' Ideal 
linear ' ,  ' LS  weights ' ) ; 


GdBavg=GdBavg+GdB ; 

GdBerravg=GdBerravg+GdBerr; 

GdBerrlinavg=GdBerrlinavg+GdBerrlin; 

GdBref avg=GdBref avg+GdBref ; 

GdBlsavg=GdBlsavg+GdBls ; 

GdBLSavg=GdBLSavg+GdBLS ; 

GdBLSlavg=GdBLSlavg+GdBLSl ; 

GdBLS2avg=GdBLS2avg+GdBLS2 ; 

GdBLS3avg=GdBLS3avg+GdBLS3 ; 

GdBLS  4 avg=GdBLS4 avg+GdBLS  4 ; 

GdBLS5avg=GdBLS5avg+GdBLS5 ; 

end  %%%%%  End  of  loop  for  iterated  computations  (average  Gain) 

o  o  o  o  o  o  o 
o  o  o  o  o  o  o 


GdBavg=GdBavg/NumIter; 

GdBerravg=GdBerravg/NumIter; 

GdBerrlinavg=GdBerrlinavg/NumIter; 

GdBrefavg=GdBrefavg/NumIter; 

GdBlsavg=GdBlsavg/NumIter ; 

GdBLSavg=GdBLSavg/NumIter ; 

GdBLS 1 avg=GdBLSl avg/Numlter; 

GdBLS2avg=GdBLS2avg/NumIter; 

GdBLS3avg=GdBLS3avg/NumIter; 

GdBLS4avg=GdBLS4avg/NumIter ; 

GdBLS5avg=GdBLS5avg/NumIter ; 

figure ( 30 )  ; 
theta=-90 : 90 ; 

plot ( theta  f  GdBavg ( : , phi_ang+90  +  l ) ,  1 Linewidth ' ,  2 ) ; 
hold  on; 

plot ( theta ,  GdBerravg ( : , phi_ang+90+l ) , ' r- . ' , ' Linewidth ' ,2) ; 
hold  on; 

plot (theta, GdBerrlinavg ( : , phi_ang+90+l ) , ’ m : ’ , ’ Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBref avg ( : , phi_ang+90+l ) , 1 g- ' , 1 Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBlsavg ( : , phi_ang+90+l ) , ’ c- ' , 1 Linewidth ' , 2 ) ; 
grid  on 

legend (' Correct Wrongly  estimated Assumed  perfect  linear Ideal 
linear ' , ' LS  weights ’ ) ; 
axis ( [-85  85  -50  5] ) ; 

title ('Fig. 3  :  Average  Beampattern  for  N  linear  array  elements  and 

given  \phi ' , ' Fontsize ' , 12 ) ; 

xlabel ( ' \ theta  (degrees ) ' , ' Fontsize ' , 12 ) ; 

ylabel (' Power  Gain  (dB) , '  ,' Fontsize ', 12 ) ; 
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figure  (40)  ; 
theta=-90 : 90 ; 

plot ( theta ,  GdBavg ( : , phi_ang+90  +  l ) ,  1 Linewidth  '  ,  2 ) ; 
hold  on; 

plot ( theta , GdBerravg ( : , phi_ang+90  +  l ) ,  '  r  -  .  '  ,  ' Linewidth ' ,  2 ) ; 
hold  on; 

plot ( theta  ,  GdBref avg ( : , phi_ang+90  +  l ) ,  ' g- ' ,  ' Linewidth ' ,  2 ) ; 
hold  on; 

plot ( theta ,  GdBlsavg ( : r phi_ang+90+l ) , 1 c- ' ,  ' Linewidth ’ ,  2 ) ; 
axis ( [-85  85  -50  5] )  ; 

title ('Fig. 4  :  Average  Beampattern  for  N  linear  array 
elements ' , ' Font size ' ,12) ; 

xlabel ( ' \ theta  (degrees ) ' ,  ' Fonts ize ' , 12 ) ; 
ylabel (' Power  Gain  (dB) , '  , ' Fontsize ' , 12 ) ; 

grid  on 

legend (' Correct '  ,  ' Wrongly  estimated ' ,  ' Ideal  linear',  'LS  weights'); 

figure  ( 50 ) ; 
theta=-90 : 90 ; 

plot (theta, GdBavg ( : , phi_ang+90+l ) , ' Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBref avg ( : , phi_ang+90  +  l ) ,  ' g- '  ,  ' Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBlsavg ( : , phi_ang+90+l ) , ' c- ' , ' Linewidth ' , 2 ) ; 
hold  on; 

plot (theta, GdBLSavg ( : , phi_ang+90+l ) , ' k : ' , ' Linewidth ' , 2 ) ; 
axis ( [-85  85  -50  5]); 

title ('Fig. 5  :  Average  Beampattern  for  N  linear  array 
elements ' , ' Fontsize ' , 12 ) ; 

xlabel ( ' \ theta  (degrees )  ' ,  ' Fontsize '  ,  12 ) ; 
ylabel (' Power  Gain  (dB) , '  , ' Fontsize ', 12 ) ; 

grid  on 

legend (' Correct ',' Ideal  linear', 'LS  weights', 'LS  less  constraints') 
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Array2D.m 


o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 


Filename : 
Author : 

Date : 

Description : 


InputParameters .m 
Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  function  creates  a  GUI  for  defining  the 
characteristics  of  a  random  array 


function  answer=inputparameters ; 


prompt={ 'Give  Nx  number  of  array  elements  in  x  direction:  '  ,  .  .  . 

'Give  Ny  number  of  array  elements  in  y  direction:  ',  .  .  . 
'Generation  of  position:  (1)  :Deviation  from  perfect  linear,  (2)  : 
From  scratch' , . . . 

'Position  error  (with  respect  to  perfect  linear)  in  x  direction 
'Position  error  (with  respect  to  perfect  linear)  in  y  direction 

(%)’,... 

'Estimated  position  error  (with  respect  to  actual)  in  x 
direction  (%)  '  ,  .  .  . 

'Estimated  position  error  (with  respect  to  actual)  in  y 
direction  (%)',... 

'Elevation  angle  (theta)  (degrees) 

'Azimuth  angle  (phi)  (degrees) 

'Angle  phi  for  beampattern: ' , . . . 

'Angle  error  in  theta  (  +  -  degrees)  '  ,  .  .  . 

'Angle  error  in  phi  (+-degrees) ' , . . . 

'Number  of  iterations',}; 

name= ' Parameters  for  antenna  array'; 
numlines=l ; 

def aultanswer= { ' 10 ' , ' 1 ' , ' 1 ' , ' 20 ' , ' 20 ' , ' 0 ' , ' 0 ' , ' 30 ' , ' 45 ' , ' 45 ' , ' 0 ' , ' 0 ' , ' 2 

5’  }; 

answer=inputdlg (prompt, name, numlines, def aultanswer ) ; 

for  i=l : length (answer ) ; 

temp (i) =str2num (answer} i } ) ; 

end 

answer=temp; 
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•  Gain2D.m 


o  o  o 
o  o  o 


o  o  o 
o  o  o 


Filename : 
Author : 


Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  function  calculates  the  beampattern  gain  of 
an  array 


Gain2D . m 


o  o  o 
o  o  o 


o  o  o 
o  o  o 


Date : 

Description : 


o  o  o 
o  o  o 


o  o  o 
o  o  o 


function  Gain=gain2 (wm,xm,ym) ; 


global  b 


for  theta=-90 : 90 ; 
for  phi=-90 : 90 ; 


th=theta*pi/180 ; 
ph=phi*pi/180 ; 


F ( 90+theta+l, 90+phi  +  l) =sum (sum (conj  (wm)  . *exp ( j  *b* (xm*sin (th) *cos (ph) +ym 
*sin  (th) *sin (ph) ) ) ) ) ; 
end 

end 

Gain=abs (F) . A2 ; 
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weights2.m  : 


o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

uniform 

o  o  o 
o  o  o 


Filename : 
Author : 

Date : 

Description : 


weights2 . m 
Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  function  calculates  the  weights  for  a 
array 


function  w=weights2 (xm, ym, Theta,  Phi)  ; 


global  Im  b 

w=Im. *exp ( j  ^b^ (xm^sin (Theta) *cos (Phi) +ym*sin (Theta) * sin (Phi) ) ) ; 
weights 
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•  rand  inter  dist.m  : 


'o'6'o  Filename  : 

%%%  Author: 

o  o  o 
o  o  o 

%%%  Date: 

%%%  Description: 

deviations 


distribution 


rand_inter_dist . m 
Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  function  creates  a  random  array,  the 
from  the  ideal  array  follow  a  uniform 


function  [xx, yy] =rand_inter_dist (Nx, Ny) ; 
global  1 

Ny=Ny+l;  %  First  line  will  be  ignored 
Nx=Nx+l ; 

xx=zeros (Ny,Nx) ; 
yy=zeros (Ny,Nx) ; 

for  i=l : Ny ; 

for  j  =  1 : Nx; 
if  j— 1 

xx (i, j ) =0; 

else 

xx (i, j ) =xx (iA  j— 1)  +  (rand* 1/2  +  1/ 4 ) ; 

end 

if  i==l ; 

yy  (I#-  j )  =0; 

else 

YY (if  j ) =YY (i-1/ j )  +  (rand* 1/2  +  1/ 4)  ; 

end 

end 

end 

if  Ny~=l 

xxl=xx (2 :Nyf 2 :Nx) ; 
yyl=YY (2 :Nyf  2 :Nx)  ; 

end 

xxl=reshape (xxl,  (Nx-1) * (Ny— 1) , 1) ;  %  ignore  first  line 

yyl=reshape (yyl, (Nx-1) * (Ny-1) , 1) ; 

xx=xxl ; 
yy=yyl ; 
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CostAnalysis.m 


o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 


Filename : 
Author : 

Date : 

Description : 


CostAnalysis .m 
Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  file  is  used  for  the  calculation  of  the 
processing  and  communication  costs  and  the 
power  consumption  for  the 

centralized, the  distributed  QR  decomposition 
and  the  iterative  approach. 


clear  all 
close  all 
clc 


N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

%%  Reference 
%%  Computational  Cost 
Comp_ref=2  *na2  * (M-N/3 ) +M*N+NA2 ; 
Com  refl=2*N; 


ooooooooooooooo 


Distributed  QR  decomposition 


oooooooooooooooo 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


o  o 
o  o 

COMPUTATIONAL  COST 

o  o 
o  o 

Number  of  Operations) 

o  o 
o  o 

Note  that  there  is 

no  distinguish  between  the  type  of  operations 

o  o 
o  o 

like  additions  or 

multiplications 

o 

o 

A  matrix  (M  x  N) 

o 

o 

b  vector  (M  x  1) 

o 

o 

QR  decomposition  : 

Cqr=2*NA2* (M-N/3) 

o 

o 

Update  : 

Cu=M*N 

o 

o 

Back  Substitution: 

Cb=NA2 

Cqr=2  *na2  * (M-N/3)  ; 

Cu=M^N; 

Cb=NA2 ; 

Cub=Cu+Cb; 

Compl=Cqr+Cub; 

figure ( 1 ) ; 

semi logy (M, Compl , ’ bo- ’ , M, Cqr , ’ rp- . ' , M, Cub, ’ md-- ' , ' Linewidth ',1.5) ; 
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title (' Fig. 1  :  Computational  Cost  as  a  function  of  M  , 

(N=l 0 ) ' , 'Fontsize' , 12) ; 

xlabel (' Number  of  angles  in  beampattern  approximation ' ,  ' Fontsize ! , 12 ) ; 
ylabel (' Number  of  instructions','Fontsize',12); 

legend (' Total ',' QR  decomposition ',' Update  and  Back  Substitution'); 
grid  on; 


M=2  0 ; 

N=5 : 2  0 ; 

Cqr=2*N. A2 . * (M-N/3)  ; 

Cu=M*N; 

Cb=N . A2 ; 

Cub=Cu+Cb; 

Comp2=Cqr+Cub; 

figure ( 2 ) ; 

semi logy (N, Comp2 ,  'bo-' ,N, Cqr , ' rp- . ',N, Cub , 1 md-- ' , ' Linewidth ’,1.5) ; 
title (' Fig. 1  :  Computational  Cost  as  a  function  of  M  , 

(N=10 ) ' /  ' Fontsize ' , 12 ) ; 

xlabel (' Number  of  sensors' , 'Fontsize' , 12) ; 
ylabel (' Number  of  instructions','Fontsize',12); 

legend (’ Total ' ,  ’ QR  decomposition ' ,  ' Update  and  Back  Substitution'); 
grid  on; 


%%  COMMUNICATION  COST 

%%  Defined  as  the  number  of  data  values  we  need  to  send  (broadcasting) 
%%  Not  the  number  of  packets 

%  1st  pass  :  From  1st  sensor  to  the  last  (QR  decomposition)  : 

%  Cqr= (M+4-N/2) * (N-l) 

%  2nd  pass  :  Back  substitution  phase 
%  Cbs=2  * (N-l ) 

N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

Comm_qr= (M+2-N/2) * (N-l)  ; 

Comm_bs=2* (N-l)  ; 

Coml=Comm_qr+Comm_bs ; 

figure  ( 10 )  ; 

semilogy (M, Coml ,  ' bo- ' , M, Comm_qr r  ' rp- 

. ' , M, ones ( length (M) , 1 ) *Comm_bs ,  ' md : ' , M, ones ( length (M) , 1 ) *Com_ref 1 , 1 gs- 
. ' ,  ' Linewidth ' , 1 . 5 ) ; 

title ('Fig. 3  : Communication  Cost  as  a  function  of  M  , 

(N=10 ) ' ,  ' Fontsize ',12) ; 

xlabel (' Number  of  angles  in  beampattern  approximation ' ,  ' Fontsize ' , 12 ) ; 
ylabel (' Number  of  elements  transmitted ' ,  ' Fontsize ' , 12 ) ; 
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legend (' Total '  ,  ' QR  decomposition  ' ,  ' Update  and  Back 
Substitution' ,  'Centralized  approach' ) ; 
grid  on; 


M=2  0 ; 

N=5 : 2  0 ; 

Com_ref 2=2  *N; 

Comm_qr= (M+2-N/2) . * (N-l) ; 
Comm_bs=2* (N-l) ; 

Com2=Comm_qr+Comm_bs ; 


figure (11) ; 

semilogy (N, Com2 ,  ' bo- ' , N, Comm _ qr ,  ' rp- .  ' , N, Comm_bs ,  ' md :  ' , Com_ref 2 ,  ' gs 

. ' ,  ' Linewidth ',1.5) ; 

title ('Fig. 4  : Communication  Cost  as  a  function  of  N  , 

(M=2  0 )  ' ,  'Font size' , 12) ; 

xlabel (' Number  of  sensors ',' Fontsize ', 12 ) ; 

ylabel (' Number  of  elements  transmitted','Fontsize',12); 

legend (' Total ',' QR  decomposition  ', 'Update  and  Back 

Substitution' , 'Centralized  approach' ) ; 

grid  on; 


%%%  Power  analysis 
Ntp=200 ; 

Pi=l; 

Ptb=Ntp^Pi; 

B=32; 


%  total  Power 
%  vs  M 

Pl=Compl*Pi+Coml*B*Ptb; 

Pref l=Compl^Pi+Com_ref l^B^Ptb; 

N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

figure (20 ) ; 

semilogy (M, PI , 'bo-',M, Pref 1 , ' rp : ' , ' Linewidth ',1.5) ; 
title ( ' Fig . 2 0  : Power  analysis  as  a  function  of  M  , 

(N=l 0 ) ' , 'Fontsize' , 12) ; 

xlabel (' Number  of  angles  in  beampattern  approximation ',' Fontsize ', 12 ) 

ylabel (' Power  (number  of  Pi )',' Fontsize ', 12 ) ; 

legend (' Distributed  approach ',' Centralized  approach'); 

grid  on; 

%  vs  N 

P2=Comp2  *Pi+Com2  *B*Ptb; 
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Pref 2=Comp2  *Pi+Com_ref 2  *B*Ptb; 


M=2  0 ; 

N=5 : 2  0 ; 

figure (25) ; 

semi logy (N,  P2 ,  ' bo- ' , N, Pref 2 ,  ' rp :  !  ,  ' Linewidth ',1.5) ; 
title ( 1  Fig . 25  : Power  analysis  as  a  function  of  N  , 

(M=2  0 )  ' ,  'Font size' , 12) ; 

xlabel (' Number  of  sensors ' ,  ' Fontsize ', 12 ) ; 

ylabel (' Power  (number  of  Pi) '  , 'Fontsize', 12) ; 

legend (' Distributed  approach ' ,  ' Centralized  approach'); 

grid  on; 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-  T a -r  -h ^  h  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  2- 

ooooooooooooooo  1  LcIdLlvc  ripp  (J  cl  UI1  ooooooooooooooooooo 


N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

K=5 ; 

CP1=N* (2* (M-l/3) +K* (3*M+1) ) ; 
figure  ( 30 ) ; 

semi logy (M, CPI ,  ' bo- ' , M, Compl ,  ' rp :  '  ,  ' Linewidth ',1.5) ; 
title (' Fig . 30  : Computational  Cost  as  a  function  of  M  , 

(N=10,K=5) ' , 'Fontsize' ,12); 

xlabel (' Number  of  angles  in  beampattern  approximation ',' Fontsize ', 12 ) 
ylabel (' Number  of  instructions ',' Fontsize ', 12 ) ; 
legend (' Distributed  approach ',' Centralized  approach'); 
grid  on; 


M=2  0 ; 

N=5 : 2  0 ; 

K=4 ; 

CP2=N* (2* (M-l/3) +K* (3*M+1)  )  ; 
figure  ( 35 ) ; 

semi logy (N, CP2 , ' bo- ' , N, Comp2 , ' rp : ' , ' Linewidth ',1.5) ; 
title (' Fig . 35  : Computational  Cost  as  a  function  of  N  , 
(M=2 0 ) ' , 'Fontsize' , 12) ; 

xlabel (' Number  of  sensors' , 'Fontsize' , 12) ; 
ylabel (' Number  of  instructions ',' Fontsize ', 12 ) ; 
legend (' Distributed  approach ',' Centralized  approach'); 
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grid  on; 


M=2  0 ; 

N=10 ; 

K=1 : 8; 

CP3=N* (2* (M-l/3) +K* (3*M+1) ) ; 

Cref=2*NA2* (M-N/3 ) +M*N+NA2 ; 

figure  ( 40 )  ; 

semi logy (K, CP3 ,  ' bo- ' , K,  Cref^ones ( length  (K)  ,  1 ) ,  ' rp : ' ,  ' Linewidth ',1.5) 
title (* Fig . 40  : Computational  Cost  as  a  function  of  iterations  K  , 
(N=10,M=20) ' , ’Font size' , 12) ; 

xlabel (' Number  of  iterations ',' Fontsize ' , 12 ) ; 
ylabel (’ Number  of  instructions ’ ,  ' Fontsize 1 , 12 ) ; 
legend (’ Distributed  approach ’ ,  ' Centralized  approach'); 
grid  on; 


N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

K=5 ; 

ComParl= (K+l) *N; 
figure  ( 50 ) 


M=2  0 ; 

N=5 : 2  0 ; 

K=5 ; 

ComPar2= (K+l) *N; 

semi logy (N, ComPar2 ,  ' bo- ' , N, Com_ref 2 ,  ' rp :  ' ,  1 Linewidth  *,1.5) ; 
title (' Fig . 50  : Communication  Cost  as  a  function  of  N  , 
(M=20,K=5) ' , 'Fontsize' , 12) ; 

xlabel (' Number  of  sensors ' ,  ' Fontsize ' , 12 ) ; 
ylabel (' Number  of  elements  transmitted ' ,  ' Fontsize ' , 12 ) ; 
legend (' Distributed  approach ',' Centralized  approach'); 
grid  on; 


figure ( 60 ) ; 

M=2  0 ; 

N=10 ; 

K=1 : 8; 

ComPar3= (K+l) *N; 
Com  ref3=2*N; 
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semilogy (K, ComPar3 ,  'bo- 

'  ,  K, ones ( length (K) , 1 ) *Com_ref 3 ,  ' rp : ' ,  ' Linewidth !,1.5) ; 
title (' Fig . 60  : Communication  Cost  as  a  function  of  K  , 
(N=10,M=20)  ' ,  ' Font size  1 , 12) ; 

xlabel ( 1  Number  of  iterations ' ,  ' Fontsize ', 12 )  ; 
ylabel (' Number  of  elements  transmitted ' ,  ' Fontsize ' , 12 ) ; 
legend (' Distributed  approach ' ,  ' Centralized  approach'); 
grid  on; 


%%%  Power  Analysis 
Ntp=200 ; 

Pi=l; 

Ptb=Ntp*Pi ; 

B=32 ; 


N=10;  %  Number  of  sensors 

M=10:30;  %  Number  of  angles 

K=4 ; 

Parl=CPl*Pi+ComParl*B*Ptb; 

Pref l=Compl*Pi+Com_ref l*B*Ptb; 

figure ( 70 )  ; 

plot (M, Pari,  ' bo- ' , M, Pref 1 , ' rp : ' ,  ' Linewidth ’,1.5) ; 
title ( ' Fig . 7 0  : Power  analysis  as  a  function  of  M  , 

(N=10 ) ’ , ' Fontsize ' , 12 ) ; 

xlabel (’ Number  of  angles  in  beampattern  approximation ',' Fontsize ', 12 ) 

ylabel (' Power  (number  of  Pi )',' Fontsize ', 12 )  ; 

legend (’ Distributed  approach ',' Centralized  approach'); 

grid  on; 


M=2  0 ; 

N=5 : 2  0 ; 

K=4 ; 

Par2=CP2  *Pi+ComPar2  *B*Ptb; 

Pref 2=Comp2  *Pi+Com_ref 2  *B*Ptb; 

figure ( 80 ) ; 

semilogy (N, Par 2 , ' bo- ' , N, Pref 2 , ' rp : ' , ' Linewidth ’,1.5) ; 
title (' Fig . 80  : Power  analysis  as  a  function  of  N  , 

(M=2  0 )  ' ,  'Fontsize' ,  12)  ; 

xlabel (' Number  of  sensors ',' Fontsize ', 12 ) ; 

ylabel (' Power  (number  of  Pi )',' Fontsize ', 12 ) ; 

legend (' Distributed  approach ',' Centralized  approach'); 

grid  on; 


M=2  0 ; 
N=1 0 ; 
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K=1 : 8; 


Par3=CP3*Pi+ComPar3*B*Ptb; 

Cqr=2*N. A2 . * (M-N/3) ; 

Cu=M*N; 

Cb=N . A2 ; 

Cub=Cu+Cb; 

Comp3=Cqr+Cub ; 

Pref 3=Comp3*Pi+Com_ref 3*B*Ptb; 
figure  (  90 ) ; 

plot(K,Par3,  ' bo- ' , K,  ones ( length (K) , 1 ) *Pref 3 ,  '  rp  :  '  ,  ' Linewidth ' r  1 . 5 ) 
title (' Fig . 90  : Power  analysis  as  a  function  of  K  , 

(M=2 0 ) ' ,  'Font size' , 12) ; 

xlabel (' Number  of  iterations ' , ' Fontsize ' , 12 ) ; 

ylabel (' Power  (number  of  Pi )',' Fontsize ', 12 ) ; 

legend (' Distributed  approach ' ,  ' Centralized  approach'); 

grid  on; 
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•  Iterative.m 


o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 

o  o  o 
o  o  o 


Filename : 
Author : 

Date : 

Description : 


Iterative . m 
Nikolaos  Papalexidis 
Hellenic  Air  Force 
June  2007 

This  file  is  used  for  the  implementation  of  the 
distributed  iterative  solution  of  the  LS  problem 


xpos=x; 

ypos=y; 

A=D5 1 ; 
q=FDES5 ' ; 

[Mf  N] =size  (A)  ; 

clear  zl  z2  z3  z4  z5  z6  z7  z8  zz  n  e 

21-f]  ; 

Z2=[]  ; 

Z3=[]  ; 

Z4=[]  ; 

Z5=[]  ; 

Z6=[]  ; 

Z7-[]; 

Z8=[]  ; 

z  1  =  0  ; 
z  2  =  0  ; 
z3=0 ; 
z  4  =  0  ; 
z5=0 ; 
z6=0; 
z  7  =  0  ; 
z  8  =  0  ; 

A1=A(:,1)  ; 

A2=A  2)  ; 

A3=A(:f  3)  ; 

A4=A ( :,4); 

A5=A  (:  ,  5 )  ; 

A6=A ( : , 6)  ; 

A7=A(:f  7)  ; 

A8=A(:f  8)  ; 

zzl=zeros (N, 1) ; 
zzl (1) =z 1 / 
zz2=zeros (N, 1) ; 
zz2 (2 ) =z2 ; 
zz3=zeros (N, 1 ) ; 
zz3 (3) =z3 ; 
zz4=zeros (N, 1) ; 
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zz4 (4) =  z4; 
zz5=zeros (N, 1 ) ; 
zz5 (5) =z5; 
zz  6=zeros (N, 1 )  ; 
zz  6 (6) =  z  6 ; 
zz7=zeros (N, 1) ; 
zz7 ( 7 ) =  z  7 ; 
zz8=zeros (N, 1) ; 
z  z  8  ( 8 )  =  z  8  ; 

s  1  =  0 ; 
s  2  =  0  ; 
s3=0 ; 
s  4  =  0 ; 
s5=0  ; 
s6=0 ; 
s  7  =  0 ; 
s  8  =  0  ; 

zz=zzl+zz2+zz3+zz4+zz5+zz6+zz7+zz8; 
rO=A*  zz ( : , 1 ) -q; 

Niter=6;  %  Number  of  iterations 
r=zeros (M, Niter*N) ;  %  residual 

r ( : , 1) =r0 ; 
rhat=rO ; 
z=A\q; 


for  k=l : Niter ; 
m= ( k— 1 ) *8  +  1; 


sl=-Al\rhat ; 
rhat=rhat+Al*sl ; 
r ( : f m+1 ) =rhat ; 
zl=zl+sl; 

Zl=  [Z1  zl] ; 

s2=-A2 \rhat ; 
rhat=rhat+A2  *s2 ; 
r ( : , m+2 ) =rhat ; 
z2=z2+s2; 

Z2  =  [ Z2  z2 ]  ; 

s3=-A3\rhat ; 
rhat=rhat+A3*s3 ; 
r ( : , m+3) =rhat ; 
z3=z3+s3 ; 

Z3= [ Z3  z3] ; 

s4=-A4\rhat ; 
rhat=rhat+A4*s4 ; 
r ( : , m+4 ) =rhat ; 
z4=z4+s4; 

Z4  =  [  Z4  z 4 ]  ; 
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s5=-A5\rhat ; 
rhat=rhat+A5*s5 ; 
r ( : , m+5) =rhat ; 
z5=z5+s5 ; 

Z5=  [  Z5  z5 ] ; ; 

s6=-A6\rhat; 
rhat=rhat+A6*s  6 ; 
r ( : , m+6 ) =rhat ; 
z  6=z  6  +  s  6 ; 

Z  6= [ Z  6  z  6 ] ; 

s7=-A7\rhat ; 
rhat=rhat+A7*s7 ; 
r  (  :  , m+7 ) =rhat ; 
z7=z7+s7; 

Z7= [ Z7  z7] ; 

s8=-A8\rhat  ; 
rhat=rhat+A8^s8 ; 
r ( : , m+8 ) =rhat  ; 
z8=z8+s8; 

Z8= [Z8  z8]  ; 

zz ( :  ,  k)  =  [zl; z2; z3; z4; z5; z6; z7; z8]  ; 

end 

for  i=l : Niter^N+1 ; 

n ( i ) =norm ( r ( : , i ) ) ; 

end 

R=norm (A* z-q) ; 

figure  ( 1 ) ; 
tl=l : length  ( Zl ) ; 

Nl=length  ( tl ) ; 

subplot (3,3,1)  ; 

plot ( tl , ones (N1 , 1 ) *real  ( z  ( 1 ) ) ,  1 r- ' , tl , real  ( Zl ) ,  1  b-  .  ’  )  ; 
title ('Real  of  w_l ' , 1 Fontsize ' , 12 ) ; 
grid  on; 

subplot (3,3,2) ; 

plot ( tl , ones (N1 , 1 ) *real ( z (2 ) ) ,  ' r- ' , tl , real ( Z  2 ) ,  ' b- .  ' ) ; 
title ('Real  of  w_2 ',' Fontsize ', 12 ) ; 
grid  on; 

subplot (3,3,3) ; 

plot ( tl , ones (N1 , 1 ) *real ( z ( 3 ) ) ,  ' r- ' , tl , real ( Z  3 ) ,  ' b- .  ' ) ; 
title ('Real  of  w_3 ',' Fontsize ', 12 ) ; 
grid  on; 
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subplot (3,3,4) ; 

plot  ( tl , ones (N1 , 1 ) *real ( z ( 4 ) ) ,  ' r- ' , tl , real ( Z4 ) ,  ' b-  .  '  ) 
title ('Real  of  w_4 ',' Fontsize ',  12 )  ; 
grid  on; 

subplot (3,3,5) ; 

plot ( tl , ones (N1 , 1 ) *real ( z ( 5 ) ) , ' r- ' , tl , real ( Z5 ) , ' b- . ' ) 
title ('Real  of  w_5 ',' Fontsize ',  12 )  ; 
grid  on; 

subplot (3,3,6); 

plot  ( tl , ones (N1 , 1 ) *real ( z ( 6) ) ,  ' r- ' , tl , real ( Z 6 ) ,  ' b- .  ' ) 
title ('Real  of  w_6 ',' Fontsize ',  12 )  ; 
grid  on; 

subplot (3,3,7) ; 

plot ( tl , ones (N1 , 1 ) *real ( z ( 7 ) ) , ' r- ' , tl , real ( Z7 ) , ' b- . ' ) 
title ('Real  of  w_7 ',' Fontsize ', 12 ) ; 
grid  on; 

subplot (3,3,8) ; 

plot ( tl , ones (N1 , 1 ) *real ( z ( 8 ) ) , ' r- ' , tl , real ( Z8 ) , ' b- . ' ) 
title ('Real  of  w_8 ',' Fontsize ', 12 ) ; 
grid  on; 

subplot (3,3,9); 
plot  (n) ; 
hold  on; 

plot (ones (Niter*N+l , 1 ) *R,  ' r- ' )  ; 
title ( ' Residual  Norm ' , ' Fontsize ' , 12 ) ; 
grid  on; 


figure (2 ) ; 
tl=l : length ( Z 1 ) ; 

Nl=length ( tl ) ; 

subplot (3,3,1) ; 

plot  ( tl , ones (N1 , 1 ) ^imag ( z ( 1 ) ) ,  ' r- ' , tl , imag (Zl ) ,  ' b- .  '  ) 
title (' Imaginary  of  w_l ',' Fontsize ', 12 ) ; 
grid  on; 

subplot (3,3,2) ; 

plot ( tl , ones (N1 , 1 ) ^imag (z (2) ) , ' r- ' , tl , imag ( Z2 ) , ' b- . ' ) 
title (' Imaginary  of  w_2 ',' Fontsize ', 12 ) ; 
grid  on; 

subplot (3,3,3) ; 

plot ( tl , ones (N1 , 1 ) *imag ( z ( 3 ) ) , ' r- ' , tl , imag ( Z3 ) , ' b- . ' ) 
title (' Imaginary  of  w_3 ',' Fontsize ', 12 ) ; 
grid  on; 
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subplot (3,3,4) ; 

plot  ( tl , ones (N1 , 1 ) *imag ( z ( 4 ) ) ,  * r- ' , tl , imag ( Z 4 ) ,  1  b-  .  '  )  ; 
title (' Imaginary  of  w_4 1 , ' Fontsize 1 , 12 ) ; 
grid  on; 

subplot (3,3,5) ; 

plot ( tl , ones (N1 , 1 ) *imag ( z ( 5 ) ) ,  ' r- ' , tl , imag ( Z  5 ) ,  1 b- .  ' ) ; 
title (’ Imaginary  of  w_5 Fontsize ’,  12 )  ; 
grid  on; 

subplot (3,3,6); 

plot ( tl , ones (N1 , 1 ) *imag ( z ( 6) ) ,  ’r-',tl, imag ( Z  6 ) ,  1 b- .  ' ) ; 
title (' Imaginary  of  w_6 Fontsize 1 , 12 ) ; 
grid  on; 

subplot (3,3,7) ; 

plot ( tl , ones (N1 , 1 ) *imag ( z ( 7 ) ) ,  ’ r- ' , tl , imag ( Z  7 ) ,  1 b- .  '  ) ; 
title (' Imaginary  of  w_7 Fontsize ’,  12 )  ; 
grid  on; 

subplot (3,3,8) ; 

plot ( tl , ones (N1 , 1 ) *imag ( z ( 8 ) ) ,  ’ r- ' , tl , imag ( Z  8 ) ,  1 b- .  '  ) ; 
title (' Imaginary  of  w_8 Fontsize 1 , 12 ) ; 
grid  on; 

subplot (3,3,9); 
plot  (n) ; 
hold  on; 

plot (ones (Niter*N+l,l) *R,  1 r- ' )  ; 
title ( 1  Residual  Norm ’ ,  ’ Fontsize ' , 12 ) ; 
grid  on; 


figure ( 3 ) ; 
tl=l : length ( Z 1 ) ; 

Nl=length ( tl ) ; 

subplot (3,3,1) ; 

plot ( tl , ones (N1 , 1 ) *abs ( z ( 1 ) ) , 1 r- ’ , tl , abs (Zl ) , ' b- . ’) ; 
title (’ Magnitude  of  w_l Fontsize 1 , 12 ) ; 
grid  on; 

subplot (3,3,2) ; 

plot  ( tl , ones (N1 , 1 ) *abs ( z (2 ) ) ,  ,r-,,tl,abs(Z2) ,  ' b- . 1 ) ; 
title (' Magnitude  of  w_2 Fontsize ’, 12 ) ; 
grid  on; 

subplot (3,3,3) ; 

plot  ( tl , ones (N1 , 1 ) *abs ( z ( 3 ) ) ,  1 r- ’ , tl , abs (Z3 ) ,  ' b- .  ') ; 
title (’ Magnitude  of  w_3 Fontsize ’, 12 ) ; 
grid  on; 
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subplot (3,3,4) ; 

plot  ( tl , ones (N1 , 1 ) *abs ( z ( 4 ) ) ,  ,r-,,tl,abs(Z4) ,  '  b- . 1 ) ; 
title (’ Magnitude  of  w_4 ' , ' Fontsize 1 , 12 ) ; 
grid  on; 

subplot (3,3,5) ; 

plot  ( tl , ones (N1 , 1 ) *abs ( z ( 5 ) ) ,  1 r- ' , tl , abs (Z5 ) ,  '  b- .  '  ) ; 
title (' Magnitude  of  w_5 Fontsize ',  12 )  ; 
grid  on; 

subplot (3,3,6); 

plot  ( tl , ones (N1 , 1 ) *abs ( z ( 6) ) ,  ’r-',tl,abs(Z6) ,  1 b- . 1 ) ; 
title (' Magnitude  of  w_6 Fontsize ',  12 )  ; 
grid  on; 

subplot (3,3,7) ; 

plot ( tl , ones (N1 , 1 ) *abs ( z ( 7 ) ) , ,r-,,tl,abs(Z7) , 1 b- . 1 ) ; 
title (' Magnitude  of  w_7 Fontsize ’, 12 ) ; 
grid  on; 

subplot (3,3,8) ; 

plot  ( tl , ones (N1 , 1 ) *abs ( z ( 8 ) ) ,  'r-',tl,abs(Z8) ,  f  b- . ’) ; 
title (' Magnitude  of  w_8 Fontsize ! , 12 ) ; 
grid  on; 

subplot (3,3,9); 
plot  (n) ; 
hold  on; 

plot (ones (Niter*N+l,l) *R,  1 r- ' )  ; 
title ( ' Residual  Norm ’ , ’ Fontsize ' , 12 ) ; 
grid  on; 


figure ( 4 ) ; 

semi logy (n, 'bo- ' ) ; 

hold  on; 

semi logy (ones (Niter*N+l,l) *R, 1 rp- 1 ) ; 

title ( ' Residual  Norm' , ' Fontsize ' , 12 ) ; 

xlabel ( 1  Number  of  local  iterations Fontsize ', 12 )  ; 

ylabel ( 1  Residual  Norm' ,  ' Fontsize 1 , 12 )  ; 

grid  on; 


figure ( 5 ) ; 
for  i=l : Niter ; 

z t emp= [ Z 1 ( 1 ) ; Z 2 ( i ) ; Z 3 ( i ) ; Z 4 ( i ) ; Z 5 ( i ) ; Z 6 ( i ) ;Z7 (i) ; Z 8 ( i )  ] 
e  ( i ) =norm ( z-ztemp) ; 

end 

figure ( 6) 

semi logy (e, ' bo- ' ) ; 

title ('Norm  of  the  weight  error Fontsize ', 12 ) ; 
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xlabel (' Number  of  complete  iterations 1 ,  1 Fontsize 1 , 12 ) ; 
ylabel('Norm  of  the  error  between  approximate  and  actual 
weights  ' ,  1 Fontsize ' , 12 ) ; 
grid  on; 


no 
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